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.process.FloatProcessor; 012import imagingbook.common.geometry.basic.Pnt2d; 013import imagingbook.common.geometry.basic.Pnt2d.PntDouble; 014import imagingbook.common.geometry.mappings.linear.ProjectiveMapping2D; 015import imagingbook.common.ij.IjUtils; 016import imagingbook.common.math.Matrix; 017 018/** 019 * Lucas-Kanade (forward-additive) matcher [1], as described in [2]. Also called the "forward-additive" algorithm. See 020 * Sec. 24.2 (Alg. 24.1) of [2] for additional details. This version assumes that the origin of R is at its center! 021 * <p> 022 * [1] B. D. Lucas and T. Kanade. "An iterative image registration technique with an application to stereo vision". In 023 * Proceedings of the 7th International Joint Conference on Artificial Intelligence IJCAI’81, pp. 674–679, Vancouver, BC 024 * (1981). 025 * <br> 026 * [2] Simon Baker and Iain Matthews, "Lucas-Kanade 20 Years On: A Unifying Framework: Part 1", CMU-RI-TR-02-16 (2002) 027 * <br> 028 * [3] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 029 * (2022). 030 * 031 * @author WB 032 * @version 2014/02/08 033 */ 034public class LucasKanadeForwardMatcher extends LucasKanadeMatcher { 035 036 private int n; // number of warp parameters 037 private FloatProcessor Ix, Iy; // gradient of the search image 038 private double qmag = Double.MAX_VALUE; // magnitude of parameter difference vector 039 private double sqrError = Double.MAX_VALUE; // squared sum of differences between I and R 040 private double[][][] S; // S[u][u][n] = the steepest descent image for dimension n at pos. u,v (same size as R) 041 042 /** 043 * Constructor. 044 * @param I the search image (of type {@link FloatProcessor}). 045 * @param R the reference image (of type {@link FloatProcessor}) 046 * @param params a parameter object of type {@link LucasKanadeMatcher.Parameters}. 047 */ 048 public LucasKanadeForwardMatcher(FloatProcessor I, FloatProcessor R, Parameters params) { 049 super(I, R, params); 050 } 051 052 /** 053 * Constructor using default parameters. 054 * @param I the search image (of type {@link FloatProcessor}). 055 * @param R the reference image (of type {@link FloatProcessor}) 056 */ 057 public LucasKanadeForwardMatcher(FloatProcessor I, FloatProcessor R) { 058 this(I, R, new Parameters()); 059 } 060 061 @Override 062 public boolean hasConverged() { 063 return (qmag < params.tolerance); 064 } 065 066 @Override 067 public double getRmsError() { 068 return Math.sqrt(sqrError); 069 } 070 071 private void initializeMatch(ProjectiveMapping2D Tinit) { 072 n = getParameters(Tinit).length; 073 Ix = gradientX(I); 074 Iy = gradientY(I); 075 iteration = 0; 076 } 077 078 @Override 079 public ProjectiveMapping2D iterateOnce(ProjectiveMapping2D Tp) { 080 if (iteration < 0) { 081 initializeMatch(Tp); 082 } 083 iteration = iteration + 1; 084 double[][] H = new double[n][n]; // n x n cumulated Hessian matrix 085 double[] dp = new double[n]; // n-dim vector \delta_p = 0 086 sqrError = 0; 087 088 if (params.showSteepestDescentImages) { 089 S = new double[wR][hR][]; // S[u][v] holds a double vector of length n 090 } 091 092 // for all positions (u,v) in R do 093 for (int u = 0; u < wR; u++) { 094 for (int v = 0; v < hR; v++) { 095 Pnt2d x = PntDouble.from(u - xc, v - yc); // position w.r.t. the center of R 096 Pnt2d xT = Tp.applyTo(x); // warp x -> x' 097 098 double gx = Ix.getInterpolatedValue(xT.getX(), xT.getY()); // interpolate the gradient at pos. xw 099 double gy = Iy.getInterpolatedValue(xT.getX(), xT.getY()); 100 101 double[] G = new double[] {gx, gy}; // interpolated gradient vector 102 103 // Step 4: Evaluate the Jacobian of the warp W at position x 104 double[][] J = Tp.getJacobian(x); 105 106 // Step 5: compute the steepest descent image: 107 double[] sx = Matrix.multiply(G, J); // I_steepest = gradI(xy) * dW/dp(xy) is a n-dim vector 108 109 if (params.showSteepestDescentImages && iteration == 1) { 110 S[u][v] = sx; 111 } 112 113 // Step 6: Update the Hessian matrix 114 double[][] Hx = Matrix.outerProduct(sx, sx); 115 H = Matrix.add(H, Hx); 116 117 // Step 7: Compute sum_x [gradI*dW/dp]^T (R(x)-I(W(x;p))] = Isteepest^T * Exy] 118 double d = R.getf(u, v) - I.getInterpolatedValue(xT.getX(), xT.getY()); 119 sqrError = sqrError + d * d; 120 121 dp = Matrix.add(dp, Matrix.multiply(d, sx)); 122 } 123 } 124 125 if (params.showHessians && iteration == 1) { 126 IjUtils.createImage("H", H).show(); 127 IjUtils.createImage("Hi", Matrix.inverse(H)).show(); 128 } 129 130 // Step 8: Compute delta_p using Equation 10 131 // double[][] Hi = Matrix.invert(H); 132 // if (Hi == null) { 133 // IJ.log("could not invert Hessian matrix!"); 134 // return null; 135 // } 136 // Step 9: update the parameter vector 137 // double[] dp = Matrix.multiply(Hi, S); 138 139 // Step 8/9 (alternative) 140 double[] qopt = Matrix.solve(H, dp); 141 if (qopt == null) { // this should not happen 142 throw new RuntimeException(this.getClass().getName() + ": Encountered singular Hessian matrix!"); 143 //return null; 144 } 145 146 double[] p = Matrix.add(getParameters(Tp), qopt); 147 qmag = Matrix.normL2squared(qopt); 148 149 if (params.showSteepestDescentImages && iteration == 1) { 150 showSteepestDescentImages(S); 151 } 152 153// return ProjectiveMapping2D.fromParameters(p); 154 return toProjectiveMap(p); 155 } 156}