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-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.common.image.matching.lucaskanade; 010 011import ij.ImagePlus; 012import ij.process.FloatProcessor; 013import imagingbook.common.geometry.basic.Pnt2d; 014import imagingbook.common.geometry.basic.Pnt2d.PntDouble; 015import imagingbook.common.geometry.mappings.linear.AffineMapping2D; 016import imagingbook.common.geometry.mappings.linear.ProjectiveMapping2D; 017import imagingbook.common.geometry.mappings.linear.Translation2D; 018import imagingbook.common.math.Matrix; 019import imagingbook.common.util.ParameterBundle; 020 021/** 022 * <p> 023 * This is the common super-class for different variants of the Lucas-Kanade matcher [1]. See Ch. 24 of [2] for 024 * additional details. 025 * </p> 026 * <p> 027 * [1] B. D. Lucas and T. Kanade. "An iterative image registration technique with an application to stereo vision". In 028 * Proceedings of the 7th International Joint Conference on Artificial Intelligence IJCAI’81, pp. 674–679, Vancouver, BC 029 * (1981). 030 * <br> 031 * [2] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 032 * (2022). 033 * </p> 034 * 035 * @author WB 036 * @version 2022/12/16 037 */ 038public abstract class LucasKanadeMatcher { 039 040 /** 041 * Default parameters for the containing class and its sub-classes; a (usually modified) instance of this class is 042 * passed to the constructor of a non-abstract sub-class. 043 */ 044 public static class Parameters implements ParameterBundle<LucasKanadeMatcher> { 045 /** Convergence limit */ 046 public double tolerance = 0.00001; 047 /** Maximum number of iterations */ 048 public int maxIterations = 100; 049 /** Set true to output debug information */ 050 public boolean debug = false; 051 /** Set true to display the steepest-descent images */ 052 public boolean showSteepestDescentImages = false; 053 /** Set true to display the Hessian matrices */ 054 public boolean showHessians = false; 055 } 056 057 final FloatProcessor I; // search image 058 final FloatProcessor R; // reference image 059 final Parameters params; // parameter object 060 061 final int wR, hR; // width/height of R 062 final double xc, yc; // center (origin) of R 063 064 int iteration = -1; 065 066 /** 067 * Constructor. 068 * @param I the search image (of type {@link FloatProcessor}) 069 * @param R the reference image (of type {@link FloatProcessor}) 070 * @param params a parameter object (of type {@link LucasKanadeMatcher.Parameters}) 071 */ 072 LucasKanadeMatcher(FloatProcessor I, FloatProcessor R, Parameters params) { 073 this.I = I; // search image 074 this.R = R; // reference image 075 this.params = params; 076 wR = R.getWidth(); 077 hR = R.getHeight(); 078 xc = 0.5 * (wR - 1); 079 yc = 0.5 * (hR - 1); 080 } 081 082 /** 083 * Calculates the projective transformation that maps the reference image R (centered at the origin) to some other 084 * quad Q. 085 * 086 * @param Q an arbitrary quad (should be inside the search image I) 087 * @return the transformation from R's bounding rectangle to Q 088 */ 089 public ProjectiveMapping2D getReferenceMappingTo(Pnt2d[] Q) { // TODO: move to plugin where used 090 Pnt2d[] Rpts = this.getReferencePoints(); 091 return ProjectiveMapping2D.fromPoints(Rpts, Q); 092 } 093 094 /** 095 * Returns the corner points of the bounding rectangle of R, centered at the origin. 096 * @return the corner points of the bounding rectangle of R 097 */ 098 public Pnt2d[] getReferencePoints() { 099 double xmin = -xc; 100 double xmax = -xc + wR - 1; 101 double ymin = -yc; 102 double ymax = -yc + hR - 1; 103 Pnt2d[] pts = new Pnt2d[4]; 104 pts[0] = PntDouble.from(xmin, ymin); 105 pts[1] = PntDouble.from(xmax, ymin); 106 pts[2] = PntDouble.from(xmax, ymax); 107 pts[3] = PntDouble.from(xmin, ymax); 108 return pts; 109 } 110 111 /** 112 * Performs the full optimization on the given image pair (I, R). 113 * 114 * @param Tinit the transformation from the reference image R to the initial search patch, assuming that R is 115 * centered at the coordinate origin! 116 * @return the transformation to the best-matching patch in the search image I (again assuming that R is centered at 117 * the coordinate origin) or null if no match was found. 118 */ 119 public ProjectiveMapping2D getMatch(ProjectiveMapping2D Tinit) { 120 ProjectiveMapping2D Tp = Tinit; 121 do { 122 Tp = iterateOnce(Tp); // to be implemented by sub-classes 123 } while (Tp != null && !hasConverged() && getIteration() < params.maxIterations); 124 return Tp; 125 } 126 127 /** 128 * Performs a single matching iteration on the given image pair (I, R). 129 * 130 * @param Tp the warp transformation from the reference image R to the initial search patch, assuming that R is 131 * centered at the coordinate origin! 132 * @return a new warp transformation (again assuming that R is centered at the coordinate origin) or null if the 133 * iteration was unsuccessful. 134 */ 135 public abstract ProjectiveMapping2D iterateOnce(ProjectiveMapping2D Tp); 136 137 /** 138 * Checks if the matcher has converged. 139 * @return true if minimization criteria have been reached. 140 */ 141 public abstract boolean hasConverged(); 142 143 /** 144 * Measures the RMS intensity difference between the reference image R and the patch in the search image I defined 145 * by the current warp Tp. 146 * 147 * @return the RMS error under the current warp 148 */ 149 public abstract double getRmsError(); 150 151 /** 152 * Returns the current iteration number. 153 * @return the current iteration number 154 */ 155 public int getIteration() { 156 return iteration; 157 } 158 159 // ------------------------------------------------------------------------------------ 160 161 FloatProcessor gradientX(FloatProcessor fp) { 162 // Sobel-kernel for x-derivatives: 163 final float[] Hx = Matrix.multiply(1f/8, new float[] { 164 -1, 0, 1, 165 -2, 0, 2, 166 -1, 0, 1 167 }); 168 FloatProcessor fpX = (FloatProcessor) fp.duplicate(); 169 fpX.convolve(Hx, 3, 3); 170 return fpX; 171 } 172 173 FloatProcessor gradientY(FloatProcessor fp) { 174 // Sobel-kernel for y-derivatives: 175 final float[] Hy = Matrix.multiply(1f/8, new float[] { 176 -1, -2, -1, 177 0, 0, 0, 178 1, 2, 1 179 }); 180 FloatProcessor fpY = (FloatProcessor) fp.duplicate(); 181 fpY.convolve(Hy, 3, 3); 182 return fpY; 183 } 184 185 // ------------------------- utility methods -------------------------- 186 187 /* We must be precise about the corner points of a rectangle: 188 * If rectangle r = <u, v, w, h>, all integer values, then the first 189 * top-left corner point (u, v) corresponds to the center of pixel 190 * (u, v). The rectangle covers w pixels horizontally, i.e., 191 * pixel 0 = (u,v), 1 = (u+1,v), ..., w-1 = (u+w-1,v). 192 * Thus ROIs must have width/height > 1! 193 */ 194 195// @Deprecated 196// private Point[] getPoints(Rectangle2D r) { // does -1 matter? YES!!! CORRECT! 197// IJ.log("getpoints1: r = " + r.toString()); 198// double x = r.getX(); 199// double y = r.getY(); 200// double w = r.getWidth(); 201// double h = r.getHeight(); 202// Point[] pts = new Point[4]; 203// pts[0] = new Point.Double(x, y); 204// pts[1] = new Point.Double(x + w - 1, y); 205// pts[2] = new Point.Double(x + w - 1, y + h - 1); 206// pts[3] = new Point.Double(x, y + h - 1); 207// //IJ.log("getpoints1: p1-4 = " + pts[0] + ", " + pts[1] + ", " + pts[2] + ", " + pts[3]); 208// return pts; 209// } 210 211 void showSteepestDescentImages(double[][][] S) { // S[u][v][n] 212 String titlePrefix = "sd"; 213 int w = S.length; 214 int h = S[0].length; 215 int n = S[0][0].length; 216 for (int i = 0; i < n; i++) { 217 FloatProcessor sdip = new FloatProcessor(w, h); 218 for (int u = 0; u < w; u++) { 219 for (int v = 0; v < h; v++) { 220 sdip.setf(u, v, (float) S[u][v][i]); 221 } 222 } 223 (new ImagePlus(titlePrefix + i, sdip)).show(); 224 } 225 } 226 227 // ported from ProjectiveMapping2D -------------------------------- 228 229 double[] getParameters(ProjectiveMapping2D map) { 230 double[][] A = map.getTransformationMatrix(); 231 return new double[] { 232 A[0][0]-1, A[0][1], A[1][0], A[1][1]-1, A[2][0], A[2][1], A[0][2], A[1][2]}; 233 //return new double[] { a00 - 1, a01, a10, a11 - 1, a20, a21, a02, a12 }; 234 } 235 236 ProjectiveMapping2D toProjectiveMap(double[] param) { 237 if (param.length < 8) { 238 throw new IllegalArgumentException("Affine mapping requires 8 parameters"); 239 } 240 return new ProjectiveMapping2D( 241 param[0] + 1, param[1], param[6], 242 param[2], param[3] + 1, param[7], 243 param[4], param[5] ); 244 } 245 246 // --------------------- 247 248 double[] getParameters(AffineMapping2D map) { 249 double[][] A = map.getTransformationMatrix(); 250 return new double[] { 251 A[0][0] - 1, A[0][1], A[1][0], A[1][1]-1, A[0][2], A[1][2]}; 252 //return new double[] { a00 - 1, a01, a10, a11 - 1, a02, a12 }; 253 } 254 255 AffineMapping2D toAffineMap(double[] param) { 256 if (param.length < 6) { 257 throw new IllegalArgumentException("Affine mapping requires 6 parameters"); 258 } 259 return new AffineMapping2D( 260 param[0] + 1, param[1], param[4], 261 param[2], param[3] + 1, param[5]); 262 } 263 264 // -------------------- 265 266 double[] getParameters(Translation2D map) { 267 double[][] A = map.getTransformationMatrix(); 268// double[] p = new double[] {a02, a12}; 269 return new double[] {A[0][2], A[1][2]}; 270 } 271 272 Translation2D toTranslation(double[] p) { 273 if (p.length < 2) { 274 throw new IllegalArgumentException("Translation requires 2 parameters"); 275 } 276 return new Translation2D(p[0], p[1]); 277 } 278 279}