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.corners; 010 011import static imagingbook.common.math.Arithmetic.EPSILON_DOUBLE; 012import static imagingbook.common.math.Arithmetic.isZero; 013import static imagingbook.common.math.Arithmetic.sqr; 014import static imagingbook.common.math.Matrix.add; 015import static imagingbook.common.math.Matrix.distL2; 016import static imagingbook.common.math.Matrix.dotProduct; 017import static imagingbook.common.math.Matrix.multiply; 018import static imagingbook.common.math.Matrix.normL2; 019 020/** 021 * <p> 022 * Common interface for all sub-pixel maximum locator implementations. A sub-pixel maximum locator tries to interpolate 023 * the continuous position and value of a local maximum for a discrete 3x3 neighborhood. The center value in the 024 * discrete neighborhood is assumed to be a local maximum. See Sec. 6.5 and Appendix E of [1] for additional details. 025 * </p> 026 * <p> 027 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 028 * (2022). 029 * </p> 030 * 031 * @author WB 032 * @version 2022/09/11 revised instance creation 033 * @see QuadraticTaylor 034 * @see QuadraticLeastSquares 035 * @see Quartic 036 */ 037public abstract class SubpixelMaxInterpolator { 038 039 private SubpixelMaxInterpolator() {} 040 041 /** 042 * Tries to locate the sub-pixel maximum from the 9 discrete sample values (s0,...,s8) taken from a 3x3 043 * neighborhood, arranged in the following order: 044 * <pre> 045 * s4 s3 s2 046 * s5 s0 s1 047 * s6 s7 s8 048 * </pre> 049 * The center value (s0) is associated with position (0,0). 050 * 051 * @param s a vector containing 9 sample values in the order described above 052 * @return a 3-element array [x,y,z], with the estimated maximum position (x,y) and the associated max. value (z). 053 * This position is relative to the center coordinate (0,0). {@code null} is returned if the maximum position could 054 * not be located. 055 */ 056 public abstract float[] getMax(float[] s); 057 058 /** 059 * Enumeration of different {@link SubpixelMaxInterpolator} methods. 060 */ 061 public enum Method { 062 /** Second-order (quadratic) Taylor interpolation */ 063 QuadraticTaylor(SubpixelMaxInterpolator.QuadraticTaylor.getInstance()), 064 065 /** Second-order (quadratic) least-squares interpolation" */ 066 QuadraticLeastSquares(SubpixelMaxInterpolator.QuadraticLeastSquares.getInstance()), 067 068 /** Quartic interpolation */ 069 Quartic(SubpixelMaxInterpolator.Quartic.getInstance()), 070 071 /** No interpolation */ 072 None(null); 073 074 private final SubpixelMaxInterpolator inst; 075 076 private Method(SubpixelMaxInterpolator instance) { 077 this.inst = instance; 078 } 079 080 public SubpixelMaxInterpolator getInstance() { 081 return this.inst; 082 } 083 } 084 085 // ------------------------------------------------------------------------------ 086 087 /** 088 * <p> 089 * 2D interpolator using second-order Taylor expansion to find the coefficients of the quadratic interpolating 090 * polynomial 091 * </p> 092 * <pre>f(x,y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 y^2 + c_5 xy</pre> 093 * <p> 094 * The resulting function fits exactly the 5 on-axis samples (s[0], s[1], s[3], s[5], s[7]) but in general does not 095 * pass through the outer diagonal samples (s[2], s[4], s[6], s[8]). See Appendix E.2.2 (Alg. E.1, 096 * 'FindMaxQuadraticTaylor') of [1] for additional details. 097 * </p> 098 * <p> 099 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, 100 * Springer (2022). 101 * </p> 102 */ 103 public static class QuadraticTaylor extends SubpixelMaxInterpolator { 104 105 private static final QuadraticTaylor instance = new QuadraticTaylor(); 106 private QuadraticTaylor() {}; 107 public static QuadraticTaylor getInstance() { 108 return instance; 109 } 110 111 private final double[] c = new double[6]; // polynomial coefficients 112 113 @Override 114 public float[] getMax(float[] s) { 115 c[0] = s[0]; 116 c[1] = (s[1] - s[5]) / 2; 117 c[2] = (s[7] - s[3]) / 2; 118 c[3] = (s[1] - 2*s[0] + s[5]) / 2; 119 c[4] = (s[3] - 2*s[0] + s[7]) / 2; 120 c[5] = (-s[2] + s[4] - s[6] + s[8]) / 4; 121 122 double d = (4*c[3]*c[4] - sqr(c[5])); // determinant of Hessian 123 124 if (d < EPSILON_DOUBLE || c[3] >= 0) { // undefined or not a maximum 125 return null; 126 } 127 128 // max position: 129 float x = (float) ((c[2]*c[5] - 2*c[1]*c[4]) / d); 130 float y = (float) ((c[1]*c[5] - 2*c[2]*c[3]) / d); 131 // max value: 132 float z = (float) (c[0] + c[1]*x + c[2]*y + c[3]*x*x + c[4]*y*y + c[5]*x*y); 133 return new float[] {x, y, z}; 134 } 135 } 136 137 /** 138 * 2D interpolator using second-order Taylor expansion to find the coefficients 139 * of the quadratic interpolating polynomial 140 * <pre>f(x,y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 y^2 + c_5 xy</pre> 141 * The resulting function fits exactly the 5 on-axis samples 142 * (s[0], s[1], s[3], s[5], s[7]) but in general does not pass 143 * through the outer diagonal samples (s[2], s[4], s[6], s[8]). 144 * Alternate version using explicit inverse Hessian. 145 */ 146 @SuppressWarnings("unused") 147 private static class QuadraticTaylorAlt extends SubpixelMaxInterpolator { 148 149 private static final QuadraticTaylorAlt instance = new QuadraticTaylorAlt(); 150 private QuadraticTaylorAlt() {}; 151 public static QuadraticTaylorAlt getInstance() { 152 return instance; 153 } 154 155 @Override 156 public float[] getMax(float[] s) { 157 double[] g = { // the gradient 158 (s[1] - s[5]) / 2, 159 (s[7] - s[3]) / 2}; 160 161 double H00 = s[1] - 2*s[0] + s[5]; // elements of the Hessian matrix 162 double H11 = s[3] - 2*s[0] + s[7]; 163 double H01 = (s[4] + s[8] - s[2] - s[6]) / 4; 164 double d = H00 * H11 - sqr(H01); // = det(H) 165 166 if (d < EPSILON_DOUBLE || H00 >= 0) { // not a maximum (minimum or saddle point) 167 return null; 168 } 169 170 double[][] Hi = { // inverse Hessian 171 { H11 / d, -H01 / d}, 172 {-H01 / d, H00 / d}}; 173 174 double[] delta = multiply(Hi, g); 175 double[] xy = multiply(-1, delta); // maximum position 176 double q = s[0] - 0.5 * dotProduct(g, delta); // maximum value 177 return new float[] {(float) xy[0], (float) xy[1], (float) q}; 178 } 179 } 180 181 // ------------------------------------------------------------------------------ 182 183 /** 184 * <p> 185 * 2D interpolator based on least-squares fitting a quadratic polynomial 186 * </p> 187 * <pre>f(x,y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 y^2 + c_5 xy</pre> 188 * <p>to the supplied sample values. 189 * See Sec. E.2.3 (Alg. E2, 'FindMaxQuadraticLeastSquares') of [1] for 190 * additional details. 191 * </p> 192 * <p> 193 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 194 * 3rd ed, Springer (2022). 195 * </p> 196 */ 197 public static class QuadraticLeastSquares extends SubpixelMaxInterpolator { 198 199 private static final QuadraticLeastSquares instance = new QuadraticLeastSquares(); 200 private QuadraticLeastSquares() {}; 201 public static QuadraticLeastSquares getInstance() { 202 return instance; 203 } 204 205 private final double[] c = new double[6]; // polynomial coefficients 206 207 @Override 208 public float[] getMax(float[] s) { 209 c[0] = (5 * s[0] + 2 * s[1] - s[2] + 2 * s[3] - s[4] + 2 * s[5] - s[6] + 2 * s[7] - s[8]) / 9; 210 c[1] = (s[1] + s[2] - s[4] - s[5] - s[6] + s[8]) / 6; 211 c[2] = (-s[2] - s[3] - s[4] + s[6] + s[7] + s[8]) / 6; 212 c[3] = (-2 * s[0] + s[1] + s[2] - 2 * s[3] + s[4] + s[5] + s[6] - 2 * s[7] + s[8]) / 6; 213 c[4] = (-2 * s[0] - 2 * s[1] + s[2] + s[3] + s[4] - 2 * s[5] + s[6] + s[7] + s[8]) / 6; 214 c[5] = (-s[2] + s[4] - s[6] + s[8]) / 4; 215 216 double d = (4*c[3]*c[4] - sqr(c[5])); 217 218 if (d < EPSILON_DOUBLE || c[3] >= 0) { // not a maximum (minimum or saddle point) 219 return null; 220 } 221 222 // max position: 223 float x = (float) ((c[2]*c[5] - 2*c[1]*c[4]) / d); 224 float y = (float) ((c[1]*c[5] - 2*c[2]*c[3]) / d); 225 // max value: 226 float z = (float) (c[0] + c[1]*x + c[2]*y + c[3]*x*x + c[4]*y*y + c[5]*x*y); 227 return new float[] {x, y, z}; 228 } 229 } 230 231 // ------------------------------------------------------------------------------ 232 233 /** 234 * <p> 235 * 2D interpolator based on fitting a 'quartic' (i.e., 4th-order) polynomial 236 * </p> 237 * <pre> 238 * f(x,y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 y^2 + c_5 x y + c_6 x^2 y + c_7 x y^2 + c_8 x^2 y^2</pre> 239 * <p> 240 * to the supplied sample values. The interpolation function passes through 241 * all sample values. The local maximum cannot be found in closed form but 242 * is found iteratively, which is not guaranteed to succeed. 243 * See Appendix E.2.4 (Alg. E.3, 'FindMaxQuartic') of [1] for 244 * additional details. 245 * </p> 246 * <p> 247 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 248 * 3rd ed, Springer (2022). 249 * </p> 250 */ 251 public static class Quartic extends SubpixelMaxInterpolator { 252 public static final int DefaultMaxIterations = 20; // iteration limit 253 public static final double DefaulMaxDelta = 1e-6; // smallest x/y move to continue search 254 public static final double DefaultMaxRad = 1.0; // x/y search boundary (-xyLimit, +xyLimit) 255 256 private final int maxIterations; 257 private final double maxDelta; 258 private final double maxRad; 259 260 private final double[] c = new double[9]; // polynomial coefficients 261 262 private static final Quartic instance = new Quartic(); 263 264 private Quartic() { 265 this(DefaultMaxIterations, DefaulMaxDelta, DefaultMaxRad); 266 } 267 268 private Quartic(int maxIterations, double maxDelta, double maxRad) { 269 this.maxIterations = maxIterations; 270 this.maxDelta = maxDelta; 271 this.maxRad = maxRad; 272 } 273 274 public static Quartic getInstance() { 275 return instance; 276 } 277 278 public static Quartic getInstance(int maxIterations, double maxDelta, double maxRad) { 279 return new Quartic(maxIterations, maxDelta, maxRad); 280 } 281 282 // ----------------------------------------------------------- 283 284 @Override 285 public float[] getMax(float[] s) { 286 c[0] = s[0]; 287 c[1] = (s[1] - s[5]) / 2; 288 c[2] = (s[7] - s[3]) / 2; 289 c[3] = (s[1] - 2 * s[0] + s[5]) / 2; 290 c[4] = (s[3] - 2 * s[0] + s[7]) / 2; 291 c[5] = (-s[2] + s[4] - s[6] + s[8]) / 4; 292 c[6] = (-s[2] + 2 *s[3] - s[4] + s[6] - 2 * s[7] + s[8]) / 4; 293 c[7] = (-2 * s[1] + s[2] - s[4] + 2 * s[5] - s[6] + s[8]) / 4; 294 c[8] = s[0] + (-2 * s[1] + s[2] - 2 * s[3] + s[4] - 2 * s[5] + s[6] - 2 * s[7] + s[8]) / 4; 295 296 // iteratively find the max location: 297 boolean done = false; 298 int n = 0; 299 double[] xCur = {0, 0}; 300 double H00 = 0, H11 = 0, H01 = 0, d = 0; 301 while (!done && n < maxIterations && normL2(xCur) < maxRad) { 302 double[] g = this.getGradient(xCur); 303 double[][] H = this.getHessian(xCur); 304 H00 = H[0][0]; 305 H11 = H[1][1]; 306 H01 = H[0][1]; 307 d = H00 * H11 - sqr(H01); 308 if (isZero(d)) { 309// throw new RuntimeException(Quartic.class.getSimpleName() + ": zero determinant"); 310 return null; 311 } 312 double[][] Hi = { // inverse Hessian 313 { H11 / d, -H01 / d}, 314 {-H01 / d, H00 / d}}; 315 double[] xNext = add(xCur, multiply(multiply(-1, Hi), g)); 316 if (distL2(xCur, xNext) < maxDelta) { 317 done = true; 318 } 319 else { 320 xCur[0] = xNext[0]; 321 xCur[1] = xNext[1]; 322 n = n + 1; 323 } 324 } 325 326 boolean isMax = (d > 0) && (H00 < 0); // false if saddle point or relative minimum 327 328 if (done && isMax) { 329 float z = (float) getInterpolatedValue(xCur); 330 return new float[] {(float) xCur[0], (float) xCur[1], z}; 331 } 332 else { 333 return null; 334 } 335 } 336 337 private double getInterpolatedValue(double[] X) { 338 final double x = X[0]; 339 final double y = X[1]; 340 return c[0] + c[1]*x + c[2]*y + c[3]*x*x + c[4]*y*y 341 + c[5]*x*y + c[6]*x*x*y + c[7]*x*y*y + c[8]*x*x*y*y; 342 } 343 344 private double[] getGradient(double[] X) { 345 return this.getGradient(X[0], X[1]); 346 } 347 348 private double[] getGradient(double x, double y) { 349 double gx = c[1] + 2*c[3]*x + c[5]*y + 2*c[6]*x*y + c[7]*y*y + 2*c[8]*x*y*y; 350 double gy = c[2] + c[5]*x + 2*c[4]*y + 2*c[7]*x*y + c[6]*x*x + 2*c[8]*x*x*y; 351 return new double[] {gx, gy}; 352 } 353 354 private double[][] getHessian(double[] X) { 355 return this.getHessian(X[0], X[1]); 356 } 357 358 private double[][] getHessian(double x, double y) { 359 final double H00 = 2*(c[3] + y*(c[6] + c[8]*y)); 360 final double H11 = 2*(c[4] + x*(c[7] + c[8]*x)); 361 final double H01 = c[5] + 2*c[6]*x + 2*c[7]*y + 4*c[8]*x*y; 362 double[][] Hi = { 363 { H00, H01}, 364 { H01, H11}}; 365 return Hi; 366 } 367 368 } 369}