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.sift; 011 012import ij.IJ; 013import ij.process.FloatProcessor; 014import imagingbook.lib.math.Matrix; 015import imagingbook.pub.sift.scalespace.DogScaleSpace; 016import imagingbook.pub.sift.scalespace.GaussianScaleSpace; 017import imagingbook.pub.sift.scalespace.ScaleLevel; 018import imagingbook.pub.sift.scalespace.ScaleOctave; 019 020import java.util.ArrayList; 021import java.util.Collections; 022import java.util.Formatter; 023import java.util.List; 024import java.util.Locale; 025 026public class SiftDetector { 027 028 /* Default parameters; a (usually modified) instance of this class 029 * may be passed to constructor of the main class. 030 */ 031 public static class Parameters { 032 /** Set true to output debug information */ 033 public boolean DEBUG = false; 034 /** Type of neigborhood used for peak detection in 3D scale space */ 035 public NeighborhoodType nhType = NeighborhoodType.NH18; 036 /** Sampling scale (nominal smoothing level of the input image) */ 037 public double sigma_s = 0.5; 038 /** Base scale of level 0 (base smoothing) */ 039 public double sigma_0 = 1.6; 040 /** Number of octaves in Gaussian/DoG scale space */ 041 public int P = 4; 042 /** Scale steps (levels) per octave */ 043 public int Q = 3; 044 /** Min. magnitude required in DoG peak detection (abs. value) */ 045 public double t_Mag = 0.01; 046 /** Min. DoG magnitude required for extrapolated peaks (abs. value) */ 047 public double t_Peak = t_Mag; 048 /** Min. difference to all neighbors in DoG peak detection (max. 0.0005) */ 049 public double t_Extrm = 0.0; 050 /** Max. number of iterations for refining the position of a key point */ 051 public int n_Refine = 5; 052 /** Max. principal curvature ratio used to eliminate line-like structures (3..10) */ 053 public double rho_Max = 10.0; 054 /** Number of orientation bins in the feature descriptor (angular resolution) */ 055 public int n_Orient = 36; 056 /** Number of smoothing steps applied to the orientation histogram */ 057 public int n_Smooth = 2; 058 /** Min. value in orientation histogram for dominant orientations (rel. to max. entry) */ 059 public double t_DomOr = 0.8; 060 /** Spatial size factor of descriptor (relative to feature scale) */ 061 public double s_Desc = 10.0; 062 /** Number of spatial descriptor bins along each x/y axis */ 063 public int n_Spat = 4; 064 /** Number of angular descriptor bins */ 065 public int n_Angl = 8; 066 /** Max. value in normalized feature vector (0.2 recommended by Lowe) */ 067 public double t_Fclip = 0.2; 068 /** Scale factor for converting normalized features to byte values in [0,255] */ 069 public double s_Fscale = 512.0; 070 /** Set true to sort detected keypoints by response magnitude */ 071 public boolean sortKeyPoints = true; 072 } 073 074 private final Parameters params; 075 076 /** 077 * Types of 3D neighborhoods used in min/max detection 078 */ 079 public enum NeighborhoodType { 080 NH8(8), NH10(10), NH18(18), NH26(26); 081 private final int size; 082 private NeighborhoodType(int size) { 083 this.size = size; 084 } 085 } 086 087 static final float EPSILON_F = 1e-35f; 088 static final double PI2 = 2 * Math.PI; 089 090 /* non-static fields */ 091 092 private GaussianScaleSpace G; 093 private DogScaleSpace D; 094 private int nhSize; 095 096 097 /* Constructors */ 098 099 public SiftDetector(FloatProcessor fp) { 100 this(fp, new Parameters()); // uses default parameters 101 } 102 103 public SiftDetector(FloatProcessor fp, Parameters params) { 104 //normalize(fp); // was destructive, delegated to ScaleLevel.getValues() 105 this.params = params; 106 nhSize = params.nhType.size; 107 G = new GaussianScaleSpace(fp, params.sigma_s, params.sigma_0, params.P, params.Q, -1, params.Q+1); 108 D = new DogScaleSpace(G); 109 } 110 111 /** 112 * Only used for debugging: produces key points with orientation 113 * histograms attached for display purposes. 114 * 115 * @param keypoints the sequence of original key points 116 * @return a sequence of enhanced key points 117 */ 118 public List<KeyPoint> makeRichKeypoints(List<KeyPoint> keypoints) { 119 if (params.DEBUG) {IJ.log("makeSiftDescriptors...");} 120 //int cnt = 0; 121 List<KeyPoint> richKeyPoints = new ArrayList<KeyPoint>(); 122 for (KeyPoint kp : keypoints) { 123 //IJ.log(" " + (cnt++)); 124 float[] oh = getOrientationHistogram(kp); 125 smoothCircular(oh,params.n_Smooth); 126 kp.orientation_histogram = oh; // TODO: remove, for testing only!! 127 List<Double> peakOrientations = findPeakOrientationIndices(oh); 128 if (params.DEBUG && peakOrientations.size() == 0) { 129 IJ.log("insufficient orientations at " + kp.u + "/" + kp.v); 130 } 131 for (double km : peakOrientations) { 132 //for (int i=0; i<Math.min(1, peakOrientations.length); i++) { // use only 1 descriptor! 133 float phi = (float) (km * 2 * Math.PI / oh.length); // 0 <= phi < 2 PI. Should be in range +/-PI? 134 KeyPoint rkp = (KeyPoint) kp.clone(); 135 rkp.orientation = phi; 136 richKeyPoints.add(rkp); 137 } 138 } 139 if (params.DEBUG) {IJ.log("makeSiftDescriptors...done");} 140 return richKeyPoints; 141 } 142 143 /** 144 * Used for debugging/illustrations only! 145 * 146 * @param oh orientation histogram 147 * @return list of histogram indices for the peak orientations 148 */ 149 public List<Double> findPeakOrientationIndices(float[] oh) { 150 int nb = oh.length; 151 List<Double> orientIndexes = new ArrayList<Double>(nb); 152 // find the maximum entry in the orientation histogram 'oh' 153 float maxh = oh[0]; 154 for (int k = 1; k < nb; k++) { 155 if (oh[k] > maxh) 156 maxh = oh[k]; 157 } 158 159 if (maxh > 0.01f) { // ascertain minimum (non-zero) gradient energy 160 // collect all peaks > 80% of the maximum entry in 'oh' 161 float minh = maxh * (float) params.t_DomOr; 162 for (int k = 0; k < nb; k++) { // hp ~ hc ~ ha 163 // angles[k] = Float.NaN; 164 float hc = oh[k]; // center value 165 if (oh[k] > minh) { // value is min. 80% of global peak 166 float hp = oh[(k - 1 + nb) % nb]; // previous histogram 167 // value 168 float hn = oh[(k + 1) % nb]; // next histogram value 169 if (hc > hp && hc > hn) { // check if 'hc' is a local peak 170 // interpolate orientation by a quadratic function 171 // (parabola): 172 double delta = interpolateQuadratic(hp, hc, hn); 173 double k_max = (k + delta + nb) % nb; // interpolated 174 // bin index, 0 175 // <= km < nPhi 176 // double phi_max = k_max * 2 * Math.PI / nb; // 0 <= 177 // phi < 2 PI. Should be in range +/-PI? 178 orientIndexes.add(k_max); 179 } 180 } 181 } 182 } 183 return orientIndexes; 184 } 185 186 187 /** 188 * THE REAL STUFF: Creating the SIFT Descriptors 189 * 190 * @return the sequence of extracted SIFT descriptors 191 */ 192 public List<SiftDescriptor> getSiftFeatures() { 193 List<KeyPoint> keyPoints = getKeyPoints(); 194 List<SiftDescriptor> siftDescriptors = new ArrayList<SiftDescriptor>(); 195 for (KeyPoint c : keyPoints) { 196 for (double phi_d : getDominantOrientations(c)) { 197 SiftDescriptor sd = makeSiftDescriptor(c, phi_d); 198 if (sd != null) { 199 siftDescriptors.add(sd); 200 } 201 } 202 } 203 return siftDescriptors; 204 } 205 206 public List<KeyPoint> getKeyPoints() { 207 List<KeyPoint> keyPts = new ArrayList<KeyPoint>(); 208 final int P = params.P; 209 final int K = params.Q; 210 for (int p = 0; p <= P-1; p++) { // for every octave p 211 for (int q = 0; q <= K-1; q++) { // for every scale level q 212 List<KeyPoint> extrema = findExtrema(p, q); 213 for (KeyPoint e : extrema) { 214 KeyPoint c = refineKeyPosition(D, e); 215 if (c != null) { 216 keyPts.add(c); 217 } 218 } 219 } 220 } 221 if (params.sortKeyPoints) { 222 Collections.sort(keyPts); 223 } 224 return keyPts; 225 } 226 227 private List<KeyPoint> findExtrema(int p, int q) { 228 final float tMag = (float) params.t_Mag; 229 final float tExtrm = (float) params.t_Extrm; 230 ScaleOctave Dp = D.getOctave(p); 231 ScaleLevel Dpq = D.getScaleLevel(p, q); 232 int M = Dpq.getWidth(); 233 int N = Dpq.getHeight(); 234 List<KeyPoint> E = new ArrayList<KeyPoint>(); 235 float scale = (float) D.getAbsoluteScale(p, q); //D.getScaleIndexFloat(p, q); needed? 236 237 final float[][][] nh = new float[3][3][3]; // 3x3x3 neighborhood [q][u][v] 238 for (int u = 1; u <= M-2; u++) { 239 float x_real = (float) D.getRealX(p, u); // for display purposes only 240 for (int v = 1; v <= N-2; v++) { 241 float y_real = (float) D.getRealY(p, v); // for display purposes only 242 float mag = Math.abs(Dpq.getf(u, v)); 243 if (mag > tMag) { 244 Dp.getNeighborhood(q, u, v, nh); // CHANGE to use D not Dp!! 245 if (isExtremum(nh, tExtrm)) { 246 KeyPoint e = new KeyPoint(p, q, u, v, u, v, x_real, y_real, scale, mag); 247 E.add(e); 248 } 249 } 250 } 251 } 252 return E; 253 } 254 255 private KeyPoint refineKeyPosition(DogScaleSpace D, KeyPoint k) { 256 //IJ.log("++Processing " + c.toString()); 257 final int p = k.p; 258 final int q = k.q; 259 int u = k.u; 260 int v = k.v; 261 ScaleOctave Dp = D.getOctave(p); 262 //ScaleLevel Dpq = Dp.getLevel(q); 263 264 final double rhoMax = params.rho_Max; 265 266 final float[][][] nh = new float[3][3][3]; // 3x3x3 neighborhood nh[q][u][v] 267 final float[] grad = new float[3]; // gradient (dx, dy, ds) 268 final float[] d = new float[3]; // relocation displacement (dx, dy, ds) 269 final float[][] hess = new float[3][3]; // 3D Hessian matrix: 270 // | dxx dxy dxs | 271 // | dxy dyy dys | 272 // | dxs dys dss | 273 final double tPeak = params.t_Peak; 274 final int n_max = params.n_Refine; 275 276 final float aMax = (float) (sqr(rhoMax + 1) / rhoMax); 277 KeyPoint kr = null; // refined keypoint 278 boolean done = false; 279 int n = 1; 280 while (!done && n <= n_max && Dp.isInside(q, u, v)) { 281 Dp.getNeighborhood(q, u, v, nh); // TODO: CHANGE to use D instead of Dp! 282 gradient(nh, grad); // result stored in grad 283 hessian(nh, hess); // result stored in hess 284 final double detH = Matrix.determinant3x3(hess); 285 if (Math.abs(detH) < EPSILON_F) { // Hessian matrix has zero determinant? 286 done = true; // ignore this point and finish 287 } 288 else { 289 final float dxx = hess[0][0]; // extract 2x2 Hessian for calculating curvature ratio below 290 final float dxy = hess[0][1]; 291 final float dyy = hess[1][1]; 292// final float[][] invHess = Matrix.inverse3x3(hess); // deprecated 293 final float[][] invHess = Matrix.inverse(hess); // alternative (numerically stable) 294 Matrix.multiplyD(invHess, grad, d); // d <-- invHess . grad 295 Matrix.multiplyD(-1, d); // d <-- - d 296 final float xx = d[0]; // = x' 297 final float yy = d[1]; // = y' 298 299 if (Math.abs(xx) < 0.5f && Math.abs(yy) < 0.5f) { // stay at this lattice point 300 done = true; 301 final float Dpeak = nh[1][1][1] + 0.5f * (grad[0]*xx + grad[1]*yy + grad[2]*d[2]); 302 final float detHxy = dxx * dyy - dxy * dxy; 303 if (Math.abs(Dpeak) > tPeak && detHxy > 0) { // check peak magnitude 304 final float a = sqr(dxx + dyy) / detHxy; 305 if (a <= aMax) { // check curvature ratio (edgeness) 306 if (params.DEBUG) IJ.log(k.toString() + String.format(": added after %d tries, alpha = %f", n, a)); 307 kr = k; // the passed keypoint is reused for the refined keypoint, with p,q,u,v unchanged 308 kr.x = u + xx; 309 kr.y = v + yy; 310 kr.x_real = (float) D.getRealX(p, kr.x); 311 kr.y_real = (float) D.getRealY(p, kr.y); 312 kr.scale = (float) D.getAbsoluteScale(p, q); 313 } 314 } 315 } 316 else { // large displacement, move to a neighboring DoG cell at same level q by at most one unit 317 u = u + Math.min(1, Math.max(-1, Math.round(xx))); // move u by max +/-1 318 v = v + Math.min(1, Math.max(-1, Math.round(yy))); // move v by max +/-1 319 // Note: we don't move along the scale axis! 320 } 321 } 322 n = n + 1; 323 } 324 return kr; 325 } 326 327 private boolean isExtremum(float[][][] nh, float tExtrm) { 328 return isLocalMin(nh, tExtrm) || isLocalMax(nh, tExtrm); 329 } 330 331 boolean isLocalMin(final float[][][] neighborhood, float tExtrm) { 332 final float c = neighborhood[1][1][1] + tExtrm; // center value + threshold 333 if (c >= 0) return false; // local minimum must have a negative value 334 // check 8 neighbors in scale plane q 335 if (c >= neighborhood[1][0][0]) return false; 336 if (c >= neighborhood[1][1][0]) return false; 337 if (c >= neighborhood[1][2][0]) return false; 338 if (c >= neighborhood[1][0][1]) return false; 339 if (c >= neighborhood[1][2][1]) return false; 340 if (c >= neighborhood[1][0][2]) return false; 341 if (c >= neighborhood[1][1][2]) return false; 342 if (c >= neighborhood[1][2][2]) return false; 343 if (nhSize >= 10) { 344 if (c >= neighborhood[0][1][1]) return false; 345 if (c >= neighborhood[2][1][1]) return false; 346 if (nhSize >= 18) { 347 // check 4 more neighbors in scale plane q-1 348 if (c >= neighborhood[0][0][1]) return false; 349 if (c >= neighborhood[0][2][1]) return false; 350 if (c >= neighborhood[0][1][0]) return false; 351 if (c >= neighborhood[0][1][2]) return false; 352 // check 4 more neighbors in scale plane q+1 353 if (c >= neighborhood[2][0][1]) return false; 354 if (c >= neighborhood[2][2][1]) return false; 355 if (c >= neighborhood[2][1][0]) return false; 356 if (c >= neighborhood[2][1][2]) return false; 357 // optionally check 8 remaining neighbors (corners of cube): 358 if (nhSize >= 26) { 359 if (c >= neighborhood[0][0][0]) return false; 360 if (c >= neighborhood[0][2][0]) return false; 361 if (c >= neighborhood[0][0][2]) return false; 362 if (c >= neighborhood[0][2][2]) return false; 363 if (c >= neighborhood[2][0][0]) return false; 364 if (c >= neighborhood[2][2][0]) return false; 365 if (c >= neighborhood[2][0][2]) return false; 366 if (c >= neighborhood[2][2][2]) return false; 367 } 368 } 369 } 370 return true; 371 } 372 373 boolean isLocalMax(final float[][][] neighborhood, float tExtrm) { 374 final float c = neighborhood[1][1][1] - tExtrm; // center value - threshold 375 if (c <= 0) return false; // local maximum must have a positive value 376 // check 8 neighbors in scale plane q 377 if (c <= neighborhood[1][0][0]) return false; 378 if (c <= neighborhood[1][1][0]) return false; 379 if (c <= neighborhood[1][2][0]) return false; 380 if (c <= neighborhood[1][0][1]) return false; 381 if (c <= neighborhood[1][2][1]) return false; 382 if (c <= neighborhood[1][0][2]) return false; 383 if (c <= neighborhood[1][1][2]) return false; 384 if (c <= neighborhood[1][2][2]) return false; 385 if (nhSize >= 10) { 386 if (c <= neighborhood[0][1][1]) return false; 387 if (c <= neighborhood[2][1][1]) return false; 388 if (nhSize >= 18) { 389 // check 4 more neighbors in scale plane q-1 390 if (c <= neighborhood[0][0][1]) return false; 391 if (c <= neighborhood[0][2][1]) return false; 392 if (c <= neighborhood[0][1][0]) return false; 393 if (c <= neighborhood[0][1][2]) return false; 394 // check 4 more neighbors in scale plane q+1 395 if (c <= neighborhood[2][0][1]) return false; 396 if (c <= neighborhood[2][2][1]) return false; 397 if (c <= neighborhood[2][1][0]) return false; 398 if (c <= neighborhood[2][1][2]) return false; 399 // optionally check 8 remaining neighbors (corners of cube): 400 if (nhSize >= 26) { 401 if (c <= neighborhood[0][0][0]) return false; 402 if (c <= neighborhood[0][2][0]) return false; 403 if (c <= neighborhood[0][0][2]) return false; 404 if (c <= neighborhood[0][2][2]) return false; 405 if (c <= neighborhood[2][0][0]) return false; 406 if (c <= neighborhood[2][2][0]) return false; 407 if (c <= neighborhood[2][0][2]) return false; 408 if (c <= neighborhood[2][2][2]) return false; 409 } 410 } 411 } 412 return true; 413 } 414 415 /* 416 * Calculates the 3-dimensional gradient vector for the 3x3x3 neighborhood 417 * nh[s][x][y]. The result is stored in the supplied vector grad, to which a reference is 418 * returned. 419 */ 420 float[] gradient(final float[][][] nh, final float[] grad) { 421 // Note: factor 0.5f not needed, kept for clarity. 422 grad[0] = 0.5f * (nh[1][2][1] - nh[1][0][1]); // = dx 423 grad[1] = 0.5f * (nh[1][1][2] - nh[1][1][0]); // = dy 424 grad[2] = 0.5f * (nh[2][1][1] - nh[0][1][1]); // = ds 425 return grad; 426 } 427 428 /* 429 * Calculates the 3x3 Hessian matrix for the 3x3x3 neighborhood nh[s][i][j], 430 * with scale index s and spatial indices i, j. 431 * The result is stored in the supplied array hess, to which a reference is 432 * returned. 433 */ 434 float[][] hessian(final float[][][] nh, final float[][] hess) { 435 final float nh_111 = 2 * nh[1][1][1]; 436 final float dxx = nh[1][0][1] - nh_111 + nh[1][2][1]; 437 final float dyy = nh[1][1][0] - nh_111 + nh[1][1][2]; 438 final float dss = nh[0][1][1] - nh_111 + nh[2][1][1]; 439 440 final float dxy = (nh[1][2][2] - nh[1][0][2] - nh[1][2][0] + nh[1][0][0]) * 0.25f; 441 final float dxs = (nh[2][2][1] - nh[2][0][1] - nh[0][2][1] + nh[0][0][1]) * 0.25f; 442 final float dys = (nh[2][1][2] - nh[2][1][0] - nh[0][1][2] + nh[0][1][0]) * 0.25f; 443 444 hess[0][0] = dxx; hess[0][1] = dxy; hess[0][2] = dxs; 445 hess[1][0] = dxy; hess[1][1] = dyy; hess[1][2] = dys; 446 hess[2][0] = dxs; hess[2][1] = dys; hess[2][2] = dss; 447 return hess; 448 } 449 450 void printMatrix3x3(float[][] A) { 451 IJ.log(String.format(Locale.US, "{{%.6f, %.6f, %.6f},", A[0][0], A[0][1], A[0][2])); 452 IJ.log(String.format(Locale.US, " {%.6f, %.6f, %.6f},", A[1][0], A[1][1], A[1][2])); 453 IJ.log(String.format(Locale.US, " {%.6f, %.6f, %.6f}}", A[2][0], A[2][1], A[2][2])); 454 } 455 456 /* 457 * Returns a list of orientations (angles) for the keypoint c. 458 */ 459 List<Double> getDominantOrientations(KeyPoint c) { 460 float[] h_phi = getOrientationHistogram(c); 461 smoothCircular(h_phi, params.n_Smooth); 462 return findPeakOrientations(h_phi); 463 } 464 465 // Smoothes the array A in-place (i.e., destructively) in 'iterations' passes. 466 void smoothCircular(float[] X, int n_iter) { 467 final float[] H = { 0.25f, 0.5f, 0.25f }; // filter kernel 468 final int n = X.length; 469 for (int i = 1; i <= n_iter; i++) { 470 float s = X[0]; 471 float p = X[n - 1]; 472 for (int j = 0; j <= n - 2; j++) { 473 float c = X[j]; 474 X[j] = H[0] * p + H[1] * X[j] + H[2] * X[j + 1]; 475 p = c; 476 } 477 X[n - 1] = H[0] * p + H[1] * X[n - 1] + H[2] * s; 478 } 479 } 480 481 /* 482 * Extracts the peaks in the orientation histogram 'h_phi' 483 * and returns the corresponding angles in a (possibly empty) list. 484 */ 485 List<Double> findPeakOrientations(float[] h_phi) { 486 int n = h_phi.length; 487 List<Double> angles = new ArrayList<Double>(n); 488 // find the maximum entry in the orientation histogram 'h_phi' 489 float h_max = h_phi[0]; 490 for (int k = 1; k < n; k++) { 491 if (h_phi[k] > h_max) 492 h_max = h_phi[k]; 493 } 494 if (h_max > 0.01f) { // ascertain minimum (non-zero) gradient energy 495 // collect all peaks > 80% of the maximum entry in 'oh' 496 float h_min = h_max * 0.8f; 497 for (int k = 0; k < n; k++) { // hp ~ hc ~ ha 498 float hc = h_phi[k]; // center value 499 if (hc > h_min) { // value is min. 80% of global peak 500 float hp = h_phi[(k - 1 + n) % n]; // previous histogram value 501 float hn = h_phi[(k + 1) % n]; // next histogram value 502 if (hc > hp && hc > hn) { // check if 'hc' is a local peak 503 // interpolate orientation by a quadratic function (parabola): 504 double delta = interpolateQuadratic(hp, hc, hn); 505 double k_max = (k + delta + n) % n; // interpolated bin index, 0 <= km < nPhi 506 double phi_max = k_max * 2 * Math.PI / n; // 0 <= phi < 2 PI. Should be in range +/-PI? 507 angles.add(phi_max); 508 } 509 } 510 } 511 } 512 return angles; 513 } 514 515 float[] getOrientationHistogram(KeyPoint c) { 516 final int n_phi = params.n_Orient; 517 final int K = params.Q; 518 519 ScaleLevel Gpq = G.getScaleLevel(c.p, c.q); 520 final float[] h_phi = new float[n_phi]; // automatically initialized to zero 521 final int M = Gpq.getWidth(), N = Gpq.getHeight(); 522 final double x = c.x, y = c.y; 523 524 final double sigma_w = 1.5 * params.sigma_0 * Math.pow(2,(double)c.q/K); 525 final double sigma_w22 = 2 * sigma_w * sigma_w; 526 final double r_w = Math.max(1, 2.5 * sigma_w); 527 final double r_w2 = r_w * r_w; 528 529 final int u_min = Math.max((int)Math.floor(x - r_w), 1); 530 final int u_max = Math.min((int)Math.ceil(x + r_w), M-2); 531 final int v_min = Math.max((int)Math.floor(y - r_w), 1); 532 final int v_max = Math.min((int)Math.ceil(y + r_w), N-2); 533 534 double[] gradPol = new double[2]; // gradient magnitude/orientation 535 536 for (int u = u_min; u <= u_max; u++) { 537 double dx = u - x; // distance to feature's center 538 for (int v = v_min; v <= v_max; v++) { 539 double dy = v - y; 540 double r2 = dx * dx + dy * dy; // squared distance from center 541 if (r2 < r_w2) { // inside limiting circle 542 Gpq.getGradientPolar(u, v, gradPol); 543 double E = gradPol[0], phi = gradPol[1]; // gradient magnitude/orientation 544 double wG = Math.exp(-(dx*dx + dy*dy)/sigma_w22); // Gaussian weight 545 double z = E * wG; 546 double k_phi = n_phi * phi / (2 * Math.PI); // continuous histogram bin index 547 double alpha = k_phi - Math.floor(k_phi); // weight alpha 548 int k_0 = mod((int)Math.floor(k_phi), n_phi); // lower histogram index 549 int k_1 = mod(k_0 + 1, n_phi); // upper histogram index 550 h_phi[k_0] = (float) (h_phi[k_0] + (1 - alpha) * z); // distribute z into bins k_0, k_1 551 h_phi[k_1] = (float) (h_phi[k_1] + alpha * z); 552 } 553 } 554 } 555 return h_phi; 556 } 557 558 SiftDescriptor makeSiftDescriptor(KeyPoint c, double phi_d) { 559 final int p = c.p, q = c.q; 560 final double x = c.x, y = c.y; 561 final double mag = c.magnitude; 562 563 ScaleLevel Gpq = G.getScaleLevel(p, q); 564 final int M = Gpq.getWidth(), N = Gpq.getHeight(), K = G.getQ(); 565 566 double sigma_q = G.getSigma_0() * Math.pow(2, (double) q / K); // = Gpq.getAbsoluteScale() ?? decimated scale 567 double w_d = params.s_Desc * sigma_q; // descriptor width 568 double sigma_d = 0.25 * w_d; // width of Gaussian weighting function 569 double sigma_d2 = 2 * sigma_d * sigma_d; 570 double r_d = 2.5 * sigma_d; // limiting radius 571 double r_d2 = r_d * r_d; // squared limiting radius 572 573 double sc = 1.0/w_d; // scale to canonical frame 574 double sin_phi_d = Math.sin(-phi_d); // precalculated sine 575 double cos_phi_d = Math.cos(-phi_d); // precalculated cosine 576 577 // Debug(logvar(p,"p") + logvar(q,"q") + logvar(w_d,"w_d")); 578 579 int u_min = Math.max((int)Math.floor(x - r_d), 1); 580 int u_max = Math.min((int)Math.ceil(x + r_d), M - 2); 581 // Debug(logvar(x,"x") + logvar(u_min,"u_min") + logvar(u_max,"u_max")); 582 583 int v_min = Math.max((int)Math.floor(y - r_d), 1); 584 int v_max = Math.min((int)Math.ceil(y + r_d), N - 2); 585 // Debug(logvar(y,"y") + logvar(v_min,"v_min") + logvar(v_max,"v_max")); 586 587 // create the 3D orientation histogram, initialize to zero: 588 final int n_Spat = params.n_Spat; 589 final int n_Angl = params.n_Angl; 590 double[][][] h_grad = new double[n_Spat][n_Spat][n_Angl]; 591 double[] gmo = new double[2]; // gradient magnitude/orientation 592 593 for (int u = u_min; u <= u_max; u++) { 594 double dx = u - x; 595 for (int v = v_min; v <= v_max; v++) { 596 double dy = v - y; 597 double r2 = sqr(dx) + sqr(dy); // squared distance from center 598 if (r2 < r_d2) { // inside limiting circle 599 // map to the canonical coordinate frame, ii,jj \in [-1/2, +1/2]: 600 double uu = sc * (cos_phi_d * dx - sin_phi_d * dy); 601 double vv = sc * (sin_phi_d * dx + cos_phi_d * dy); 602 // calculate the gradient of Gaussian scale level (p,q) at position (u,v): 603 Gpq.getGradientPolar(u, v, gmo); 604 double E = gmo[0]; // gradient magnitude 605 double phi = gmo[1]; // gradient orientation 606 double phi_norm = mod(phi - phi_d, PI2);// normalized gradient orientation 607 double w_G = Math.exp(-r2/sigma_d2); // Gaussian weight 608 double z = E * w_G; // quantity to accumulate 609 updateGradientHistogram(h_grad, uu, vv, phi_norm, z); 610 } 611 } 612 } 613 int[] f_int = makeFeatureVector(h_grad); 614 double sigma_pq = G.getAbsoluteScale(p, q); 615 double x_real = G.getRealX(p, x); 616 double y_real = G.getRealY(p, y); 617 return new SiftDescriptor(x_real, y_real, sigma_pq, mag, phi_d, f_int); 618 } 619 620 private void updateGradientHistogram(double[][][] h_grad, double uu, double vv, double phi_norm, double z) { 621 final int n_Spat = params.n_Spat; 622 final int n_Angl = params.n_Angl; 623 624 final double ii = n_Spat * uu + 0.5 * (n_Spat - 1); // continuous spatial histogram index i' 625 final double jj = n_Spat * vv + 0.5 * (n_Spat - 1); // continuous spatial histogram index j' 626 final double kk = phi_norm * (n_Angl / PI2); // continuous orientation histogram index k' 627 628 final int i0 = (int) Math.floor(ii); 629 final int i1 = i0 + 1; 630 631 final int j0 = (int) Math.floor(jj); 632 final int j1 = j0 + 1; 633 634 final int k0 = mod((int) Math.floor(kk), n_Angl); 635 final int k1 = (k0 + 1) % n_Angl; // k0 >= 0 636 637 final double alpha0 = 1.0 - (ii - i0); 638 final double alpha1 = 1.0 - alpha0; 639 640 final double beta0 = 1.0 - (jj - j0); 641 final double beta1 = 1.0 - beta0; 642 643 final double gamma0 = 1.0 - (kk - Math.floor(kk)); 644 final double gamma1 = 1.0 - gamma0; 645 646 final int[] I = {i0, i1}; // index arrays used in loops below 647 final int[] J = {j0, j1}; 648 final int[] K = {k0, k1}; 649 650 final double[] A = {alpha0, alpha1}; 651 final double[] B = {beta0, beta1}; 652 final double[] C = {gamma0, gamma1}; 653 654 // distribute z over 8 adjacent spatial/angular bins: 655 for (int a = 0; a <= 1; a++) { 656 int i = I[a]; 657 if (i >= 0 && i < n_Spat) { 658 double wa = A[a]; 659 for (int b = 0; b <= 1; b++) { 660 int j = J[b]; 661 if (j >= 0 && j < n_Spat) { 662 double wb = B[b]; 663 for (int c = 0; c <= 1; c++) { 664 int k = K[c]; 665 double wc = C[c]; 666 h_grad[i][j][k] = h_grad[i][j][k] + z * wa * wb * wc; 667 668 } 669 } 670 } 671 } 672 } 673 } 674 675 private int[] makeFeatureVector(double[][][] h_grad) { 676 final int n_Spat = params.n_Spat; 677 final int n_Angl = params.n_Angl; 678 float[] f = new float[n_Spat * n_Spat * n_Angl]; 679 // flatten 3D histogram to a 1D vector 680 // the histogram is vectorized such that k (phi) is the fastest 681 // varying index, followed by j and i (which is the slowest index). 682 // Note: j (v) is the slowest in VLFEAT (i,j swapped) 683 int m = 0; 684 for (int i = 0; i < n_Spat; i++) { 685 for (int j = 0; j < n_Spat; j++) { 686 for (int k = 0; k < n_Angl; k++) { 687 f[m] = (float) h_grad[i][j][k]; 688 m = m + 1; 689 } 690 } 691 } 692 normalize(f); 693 clipPeaks(f, (float) params.t_Fclip); 694 normalize(f); 695 return mapToIntegers(f, (float) params.s_Fscale); 696 } 697 698 private void normalize(float[] x) { 699 final double norm = normL2(x); 700 if (norm > EPSILON_F) { 701 final float s = (float) (1.0 / norm); 702 for (int i = 0; i < x.length; i++) { 703 x[i] = s * x[i]; 704 } 705 } 706 } 707 708 private void clipPeaks(float[] x, float xmax) { 709 for (int i=0; i<x.length; i++) { 710 if (x[i] > xmax) { 711 x[i] = xmax; 712 } 713 } 714 } 715 716 private int[] mapToIntegers(float[] x, float s) { 717 int[] ivec = new int[x.length]; 718 for (int i=0; i<x.length; i++) { 719 ivec[i] = Math.round(s * x[i]); 720 } 721 return ivec; 722 } 723 724 // auxiliary methods ------------------------- 725 726// private void normalize(FloatProcessor fp) { // replaced by method in class ScaleLevel 727// double minVal = fp.getMin(); 728// double maxVal = fp.getMax(); 729// fp.add(-minVal); 730// fp.multiply(1.0 / (maxVal - minVal)); 731// } 732 733 // also in imagingbook.math.Arithmetic.java 734 private float sqr(float x) { 735 return x*x; 736 } 737 738 // also in imagingbook.math.Arithmetic.java 739 private double sqr(double x) { 740 return x*x; 741 } 742 743 // also in imagingbook.math.Arithmetic.java 744 private int mod(int a, int b) { 745 if (b == 0) 746 return a; 747 if (a * b >= 0) 748 return a - b * (a / b); 749 else 750 return a - b * (a / b - 1); 751 } 752 753 // also in imagingbook.math.Arithmetic.java 754 private double mod(double a, double n) { 755 return a - n * Math.floor(a / n); 756 } 757 758 // also in imagingbook.math.Matrix.java 759 private double normL2(float[] vec) { 760 double sum = 0; 761 for (float x : vec) { 762 sum = sum + (x * x); 763 } 764 return Math.sqrt(sum); 765 } 766 767 /* 768 * Calculate the extremal position from 3 discrete function values 769 * at x = -1, 0, +1: f(-1) = 'y1', f(0) = 'y2', f(+1) = 'y3'. 770 */ 771 private float interpolateQuadratic(float y1, float y2, float y3) { 772 float a = (y1 - 2*y2 + y3) / 2; 773 if (Math.abs(a) < 0.00001f) { 774 throw new IllegalArgumentException("quadratic interpolation failed " 775 + a + " " + y1 + " " + y2 + " " + y3); 776 } 777 float b = (y3 - y1) / 2; 778 // float c = y2; 779 float x_extrm = -b / (2*a); // extremal position 780 // float y_extrm = a*x*x + b*x + c; // extremal value (not needed) 781 return x_extrm; // x is in [-1,+1] 782 } 783 784 // ------------------------------------------- 785 786 void print(float[] arr) { 787 int linelength = 16; 788 StringBuilder sb = new StringBuilder(); 789 Formatter fm = new Formatter(sb, Locale.US); 790 for (int i = 0; i < arr.length; i++) { 791 if (i > 0 && i % linelength == 0) { 792 IJ.log(sb.toString()); 793 sb.setLength(0); 794 } 795 fm.format(" %.2f", arr[i]); 796 } 797 IJ.log(sb.toString()); 798 fm.close(); 799 } 800 801 String logvar(float x, String name) { 802 return name + " = " + String.format("%.3f ", x); 803 } 804 805 String logvar(double x, String name) { 806 return name + " = " + String.format("%.3f ", x); 807 } 808 809 String logvar(int x, String name) { 810 return name + " = " + x + " "; 811 } 812 813 void Debug(String s) { 814 IJ.log(s); 815 } 816 817 void Stop() { 818 throw new IllegalArgumentException("HALTED"); 819 } 820 821 public void printGaussianScaleSpace() { 822 G.print(); 823 } 824}