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 *******************************************************************************/
009package imagingbook.pub.lucaskanade;
010
011import ij.ImagePlus;
012import ij.process.FloatProcessor;
013import imagingbook.lib.math.Matrix;
014import imagingbook.pub.geometry.mappings.linear.ProjectiveMapping;
015
016import java.awt.geom.Point2D;
017
018/**
019 * This is the common super-class for different variants of the Lucas-Kanade
020 * matcher.
021 *  @author Wilhelm Burger
022 *  @version 2014/02/08
023 */
024public abstract class LucasKanadeMatcher {
025
026        /**
027         * Default parameters for the containing class and its sub-classes; 
028         * a (usually modified) instance of this class is passed to the constructor 
029         * of a non-abstract sub-class.
030         */
031        public static class Parameters {
032                /** Convergence limit */
033                public double tolerance = 0.00001;
034                /** Maximum number of iterations */
035                public int maxIterations = 100;
036                /** Set true to output debug information */
037                public boolean debug = false;
038                /** Set true to display the steepest-descent images */
039                public boolean showSteepestDescentImages = false;
040                /** Set true to display the Hessian matrices */
041                public boolean showHessians = false;
042        }
043        
044        final FloatProcessor I;                 // search image
045        final FloatProcessor R;                 // reference image
046        final Parameters params;                // parameter object
047        
048        final int wR, hR;                               // width/height of R
049        final double xc, yc;                    // center (origin) of R
050        
051        int iteration = -1;
052        
053        /**
054         * Creates a new Lucas-Kanade-type matcher.
055         * @param I the search image (of type {@link FloatProcessor}).
056         * @param R the reference image (of type {@link FloatProcessor})
057         * @param params a parameter object of type  {@link LucasKanadeMatcher.Parameters}.
058         */
059        protected LucasKanadeMatcher(FloatProcessor I, FloatProcessor R, Parameters params) {
060                this.I = I;     // search image
061                this.R = R;     // reference image
062                this.params = params;
063                wR = R.getWidth();
064                hR = R.getHeight();
065                xc = 0.5 * (wR - 1);
066                yc = 0.5 * (hR - 1);
067        }
068        
069        /**
070         * Calculates the transformation that maps the reference image {@code R} (centered
071         * around the origin) to the given quad Q.
072         * @param Q an arbitrary quad (should be inside the search image I);
073         * @return the transformation from {@code R}'s bounding rectangle to {@code Q}.
074         */
075        public ProjectiveMapping getReferenceMappingTo(Point2D[] Q) {
076                Point2D[] Rpts = getReferencePoints();
077                return ProjectiveMapping.makeMapping(Rpts, Q);
078        }
079        
080        /**
081         * @return the corner points of the bounding rectangle of R, centered at the origin.
082         */
083        public Point2D[] getReferencePoints() {
084                double xmin = -xc;
085                double xmax = -xc + wR - 1;
086                double ymin = -yc;
087                double ymax = -yc + hR - 1;
088                Point2D[] pts = new Point2D[4];
089                pts[0] = new Point2D.Double(xmin, ymin);
090                pts[1] = new Point2D.Double(xmax, ymin);
091                pts[2] = new Point2D.Double(xmax, ymax);
092                pts[3] = new Point2D.Double(xmin, ymax);
093                return pts;
094        }
095        
096        
097        /**
098         * Performs the full optimization on the given image pair (I, R).
099         * @param Tinit the transformation from the reference image R to
100         * the initial search patch, assuming that R is centered at the coordinate
101         * origin!
102         * @return the transformation to the best-matching patch in the search image I
103         * (again assuming that R is centered at the coordinate origin) or null if
104         * no match was found.
105         */
106        public ProjectiveMapping getMatch(ProjectiveMapping Tinit) {
107//              initializeMatch(Tinit);                 // to be implemented by sub-classes
108                ProjectiveMapping Tp = Tinit.duplicate();
109                do {
110                        Tp = iterateOnce(Tp);           // to be implemented by sub-classes
111                } while (Tp != null && !hasConverged() && getIteration() < params.maxIterations);
112                return Tp;
113        }
114        
115        
116        /**
117         * Performs a single matching iteration on the given image pair (I, R).
118         * 
119         * @param Tp the warp transformation from the reference image R to
120         * the initial search patch, assuming that R is centered at the coordinate
121         * origin! 
122         * @return a new warp transformation (again assuming that R is centered at 
123         * the coordinate origin) or null if the iteration was unsuccessful.
124         */
125        public abstract ProjectiveMapping iterateOnce(ProjectiveMapping Tp);
126        
127        
128        /**
129         * Checks if the matcher has converged.
130         * @return true if minimization criteria have been reached.
131         */
132        public abstract boolean hasConverged();
133        
134        /**
135         * Measures the RMS intensity difference between the reference image R and the 
136         * patch in the search image I defined by the current warp Tp.
137         * @return the RMS error under the current warp
138         */
139        public abstract double getRmsError();           
140        
141        public int getIteration() {
142                return iteration;
143        }
144        
145        // ------------------------------------------------------------------------------------
146        
147        protected FloatProcessor gradientX(FloatProcessor fp) {
148                // Sobel-kernel for x-derivatives:
149            final float[] Hx = Matrix.multiply(1f/8, new float[] {
150                                -1, 0, 1,
151                            -2, 0, 2,
152                            -1, 0, 1
153                            });
154            FloatProcessor fpX = (FloatProcessor) fp.duplicate();
155            fpX.convolve(Hx, 3, 3);
156            return fpX;
157        }
158        
159        protected FloatProcessor gradientY(FloatProcessor fp) {
160                // Sobel-kernel for y-derivatives:
161                final float[] Hy = Matrix.multiply(1f/8, new float[] {
162                                                -1, -2, -1,
163                                                 0,  0,  0,
164                                                 1,  2,  1
165                                                 });
166            FloatProcessor fpY = (FloatProcessor) fp.duplicate();
167            fpY.convolve(Hy, 3, 3);
168            return fpY;
169        }
170        
171        // ------------------------- utility methods --------------------------
172        
173        /* We must be precise about the corner points of a rectangle:
174         * If rectangle r = <u, v, w, h>, all integer values, then the first
175         * top-left corner point (u, v) corresponds to the center of pixel
176         * (u, v). The rectangle covers w pixels horizontally, i.e., 
177         * pixel 0 = (u,v), 1 = (u+1,v), ..., w-1 = (u+w-1,v).
178         * Thus ROIs must have width/height > 1!
179         */
180        
181//      @Deprecated
182//      private Point2D[] getPoints(Rectangle2D r) {    // does -1 matter? YES!!! CORRECT!
183//              IJ.log("getpoints1:  r = " + r.toString());
184//              double x = r.getX();
185//              double y = r.getY();
186//              double w = r.getWidth();
187//              double h = r.getHeight();
188//              Point2D[] pts = new Point2D[4];
189//              pts[0] = new Point2D.Double(x, y);
190//              pts[1] = new Point2D.Double(x + w - 1, y);
191//              pts[2] = new Point2D.Double(x + w - 1, y + h - 1);
192//              pts[3] = new Point2D.Double(x, y + h - 1);
193//              //IJ.log("getpoints1:  p1-4 = " + pts[0] + ", " + pts[1] + ", " + pts[2] + ", " + pts[3]);
194//              return pts;
195//      }
196        
197        protected void showSteepestDescentImages(double[][][] S) {      // S[u][v][n]
198                String titlePrefix = "sd";
199                int w = S.length;
200                int h = S[0].length;
201                int n = S[0][0].length;
202                for (int i = 0; i < n; i++) {
203                        FloatProcessor sdip = new FloatProcessor(w, h);
204                        for (int u = 0; u < w; u++) {
205                                for (int v = 0; v < h; v++) {
206                                        sdip.setf(u, v, (float) S[u][v][i]);
207                                }
208                        }
209                        (new ImagePlus(titlePrefix + i, sdip)).show();
210                }
211        }
212        
213
214}