Java tutorial
/* * This file is part of the LIRE project: http://lire-project.net * LIRE is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * LIRE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with LIRE; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * * We kindly ask you to refer the any or one of the following publications in * any publication mentioning or employing Lire: * * Lux Mathias, Savvas A. Chatzichristofis. Lire: Lucene Image Retrieval * An Extensible Java CBIR Library. In proceedings of the 16th ACM International * Conference on Multimedia, pp. 1085-1088, Vancouver, Canada, 2008 * URL: http://doi.acm.org/10.1145/1459359.1459577 * * Lux Mathias. Content Based Image Retrieval with LIRE. In proceedings of the * 19th ACM International Conference on Multimedia, pp. 735-738, Scottsdale, * Arizona, USA, 2011 * URL: http://dl.acm.org/citation.cfm?id=2072432 * * Mathias Lux, Oge Marques. Visual Information Retrieval using Java and LIRE * Morgan & Claypool, 2013 * URL: http://www.morganclaypool.com/doi/abs/10.2200/S00468ED1V01Y201301ICR025 * * Copyright statement: * -------------------- * (c) 2002-2013 by Mathias Lux (mathias@juggle.at) * http://www.semanticmetadata.net/lire, http://www.lire-project.net */ package net.semanticmetadata.lire.imageanalysis.features.local.sift; /** * Difference Of Gaussian detector on top of a scale space octave as described * by David Lowe \citep{Loew04}. * * BibTeX: * <pre> * @article{Lowe04, * author = {David G. Lowe}, * title = {Distinctive Image Features from Scale-Invariant Keypoints}, * journal = {International Journal of Computer Vision}, * year = {2004}, * volume = {60}, * number = {2}, * pages = {91--110}, * } * </pre> * * License: GPL * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License 2 * as published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * @author Stephan Saalfeld <saalfeld@mpi-cbg.de> * @version 0.1b */ import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.RealMatrix; import java.util.Vector; public class FloatArray2DScaleOctaveDoGDetector { /** * minimal contrast of a candidate */ //private static final float MIN_CONTRAST = 0.03f * 255.0f; //private static final float MIN_CONTRAST = 0.025f * 255.0f; private static final float MIN_CONTRAST = 0.025f; /** * maximal curvature ratio, higher values allow more edge-like responses */ private static final float MAX_CURVATURE = 10; private static final float MAX_CURVATURE_RATIO = (MAX_CURVATURE + 1) * (MAX_CURVATURE + 1) / MAX_CURVATURE; private FloatArray2DScaleOctave octave; /** * detected candidates as float triples 0=>x, 1=>y, 2=>scale index */ private Vector<float[]> candidates; public Vector<float[]> getCandidates() { return candidates; } /** * Constructor */ public FloatArray2DScaleOctaveDoGDetector() { octave = null; candidates = null; } public void run(FloatArray2DScaleOctave o) { octave = o; candidates = new Vector<float[]>(); detectCandidates(); } private void detectCandidates() { FloatArray2D[] d = octave.getD(); for (int i = d.length - 2; i >= 1; --i) { int ia = i - 1; int ib = i + 1; for (int y = d[i].height - 2; y >= 1; --y) { int r = y * d[i].width; int ra = r - d[i].width; int rb = r + d[i].width; X: for (int x = d[i].width - 2; x >= 1; --x) { int ic = i; int iac = ia; int ibc = ib; int yc = y; int rc = r; int rac = ra; int rbc = rb; int xc = x; int xa = xc - 1; int xb = xc + 1; float e111 = d[ic].data[r + xc]; // check if d(x, y, i) is an extremum // do it pipeline-friendly ;) float e000 = d[iac].data[rac + xa]; boolean isMax = e000 < e111; boolean isMin = e000 > e111; if (!(isMax || isMin)) continue; float e100 = d[iac].data[rac + xc]; isMax &= e100 < e111; isMin &= e100 > e111; if (!(isMax || isMin)) continue; float e200 = d[iac].data[rac + xb]; isMax &= e200 < e111; isMin &= e200 > e111; if (!(isMax || isMin)) continue; float e010 = d[iac].data[rc + xa]; isMax &= e010 < e111; isMin &= e010 > e111; if (!(isMax || isMin)) continue; float e110 = d[iac].data[rc + xc]; isMax &= e110 < e111; isMin &= e110 > e111; if (!(isMax || isMin)) continue; float e210 = d[iac].data[rc + xb]; isMax &= e210 < e111; isMin &= e210 > e111; if (!(isMax || isMin)) continue; float e020 = d[iac].data[rbc + xa]; isMax &= e020 < e111; isMin &= e020 > e111; if (!(isMax || isMin)) continue; float e120 = d[iac].data[rbc + xc]; isMax &= e120 < e111; isMin &= e120 > e111; if (!(isMax || isMin)) continue; float e220 = d[iac].data[rbc + xb]; isMax &= e220 < e111; isMin &= e220 > e111; if (!(isMax || isMin)) continue; float e001 = d[ic].data[rac + xa]; isMax &= e001 < e111; isMin &= e001 > e111; if (!(isMax || isMin)) continue; float e101 = d[ic].data[rac + xc]; isMax &= e101 < e111; isMin &= e101 > e111; if (!(isMax || isMin)) continue; float e201 = d[ic].data[rac + xb]; isMax &= e201 < e111; isMin &= e201 > e111; if (!(isMax || isMin)) continue; float e011 = d[ic].data[rc + xa]; isMax &= e011 < e111; isMin &= e011 > e111; if (!(isMax || isMin)) continue; float e211 = d[ic].data[rc + xb]; isMax &= e211 < e111; isMin &= e211 > e111; if (!(isMax || isMin)) continue; float e021 = d[ic].data[rbc + xa]; isMax &= e021 < e111; isMin &= e021 > e111; if (!(isMax || isMin)) continue; float e121 = d[ic].data[rbc + xc]; isMax &= e121 < e111; isMin &= e121 > e111; if (!(isMax || isMin)) continue; float e221 = d[ic].data[rbc + xb]; isMax &= e221 < e111; isMin &= e221 > e111; if (!(isMax || isMin)) continue; float e002 = d[ibc].data[rac + xa]; isMax &= e002 < e111; isMin &= e002 > e111; if (!(isMax || isMin)) continue; float e102 = d[ibc].data[rac + xc]; isMax &= e102 < e111; isMin &= e102 > e111; if (!(isMax || isMin)) continue; float e202 = d[ibc].data[rac + xb]; isMax &= e202 < e111; isMin &= e202 > e111; if (!(isMax || isMin)) continue; float e012 = d[ibc].data[rc + xa]; isMax &= e012 < e111; isMin &= e012 > e111; if (!(isMax || isMin)) continue; float e112 = d[ibc].data[rc + xc]; isMax &= e112 < e111; isMin &= e112 > e111; if (!(isMax || isMin)) continue; float e212 = d[ibc].data[rc + xb]; isMax &= e212 < e111; isMin &= e212 > e111; if (!(isMax || isMin)) continue; float e022 = d[ibc].data[rbc + xa]; isMax &= e022 < e111; isMin &= e022 > e111; if (!(isMax || isMin)) continue; float e122 = d[ibc].data[rbc + xc]; isMax &= e122 < e111; isMin &= e122 > e111; if (!(isMax || isMin)) continue; float e222 = d[ibc].data[rbc + xb]; isMax &= e222 < e111; isMin &= e222 > e111; if (!(isMax || isMin)) continue; // so it is an extremum, try to localize it with subpixel // accuracy, if it has to be moved for more than 0.5 in at // least one direction, try it again there but maximally 5 // times boolean isLocalized = false; boolean isLocalizable = true; float dx; float dy; float di; float dxx; float dyy; float dii; float dxy; float dxi; float dyi; float ox; float oy; float oi; float od = Float.MAX_VALUE; // offset square distance float fx = 0; float fy = 0; float fi = 0; int t = 5; // maximal number of re-localizations do { --t; // derive at (x, y, i) by center of difference dx = (e211 - e011) / 2.0f; dy = (e121 - e101) / 2.0f; di = (e112 - e110) / 2.0f; // create hessian at (x, y, i) by laplace float e111_2 = 2.0f * e111; dxx = e011 - e111_2 + e211; dyy = e101 - e111_2 + e121; dii = e110 - e111_2 + e112; dxy = (e221 - e021 - e201 + e001) / 4.0f; dxi = (e212 - e012 - e210 + e010) / 4.0f; dyi = (e122 - e102 - e120 + e100) / 4.0f; // invert hessian Array2DRowRealMatrix H = new Array2DRowRealMatrix( new double[][] { { (double) dxx, (double) dxy, (double) dxi }, { (double) dxy, (double) dyy, (double) dyi }, { (double) dxi, (double) dyi, (double) dii } }); RealMatrix H_inv; try { H_inv = new LUDecomposition(H).getSolver().getInverse(); } catch (RuntimeException e) { continue X; } double[][] h_inv = H_inv.getData(); // estimate the location of zero crossing being the offset of the extremum ox = -(float) h_inv[0][0] * dx - (float) h_inv[0][1] * dy - (float) h_inv[0][0] * di; oy = -(float) h_inv[1][0] * dx - (float) h_inv[1][1] * dy - (float) h_inv[1][0] * di; oi = -(float) h_inv[2][0] * dx - (float) h_inv[2][1] * dy - (float) h_inv[2][0] * di; float odc = ox * ox + oy * oy + oi * oi; if (odc < 2.0f) { if ((Math.abs(ox) > 0.5 || Math.abs(oy) > 0.5 || Math.abs(oi) > 0.5) && odc < od) { od = odc; xc = (int) Math.round((float) xc + ox); yc = (int) Math.round((float) yc + oy); ic = (int) Math.round((float) ic + oi); if (xc < 1 || yc < 1 || ic < 1 || xc > d[0].width - 2 || yc > d[0].height - 2 || ic > d.length - 2) isLocalizable = false; else { xa = xc - 1; xb = xc + 1; rc = yc * d[ic].width; rac = rc - d[ic].width; rbc = rc + d[ic].width; iac = ic - 1; ibc = ic + 1; e000 = d[iac].data[rac + xa]; e100 = d[iac].data[rac + xc]; e200 = d[iac].data[rac + xb]; e010 = d[iac].data[rc + xa]; e110 = d[iac].data[rc + xc]; e210 = d[iac].data[rc + xb]; e020 = d[iac].data[rbc + xa]; e120 = d[iac].data[rbc + xc]; e220 = d[iac].data[rbc + xb]; e001 = d[ic].data[rac + xa]; e101 = d[ic].data[rac + xc]; e201 = d[ic].data[rac + xb]; e011 = d[ic].data[rc + xa]; e111 = d[ic].data[rc + xc]; e211 = d[ic].data[rc + xb]; e021 = d[ic].data[rbc + xa]; e121 = d[ic].data[rbc + xc]; e221 = d[ic].data[rbc + xb]; e002 = d[ibc].data[rac + xa]; e102 = d[ibc].data[rac + xc]; e202 = d[ibc].data[rac + xb]; e012 = d[ibc].data[rc + xa]; e112 = d[ibc].data[rc + xc]; e212 = d[ibc].data[rc + xb]; e022 = d[ibc].data[rbc + xa]; e122 = d[ibc].data[rbc + xc]; e222 = d[ibc].data[rbc + xb]; } } else { fx = (float) xc + ox; fy = (float) yc + oy; fi = (float) ic + oi; if (fx < 0 || fy < 0 || fi < 0 || fx > d[0].width - 1 || fy > d[0].height - 1 || fi > d.length - 1) isLocalizable = false; else isLocalized = true; } } else isLocalizable = false; } while (!isLocalized && isLocalizable && t >= 0); // reject detections that could not be localized properly if (!isLocalized) { // System.err.println( "Localization failed (x: " + xc + ", y: " + yc + ", i: " + ic + ") => (ox: " + ox + ", oy: " + oy + ", oi: " + oi + ")" ); // if ( ic < 1 || ic > d.length - 2 ) // System.err.println( " Detection outside octave." ); continue; } // reject detections with very low contrast if (Math.abs(e111 + 0.5f * (dx * ox + dy * oy + di * oi)) < MIN_CONTRAST) continue; // reject edge responses float det = dxx * dyy - dxy * dxy; float trace = dxx + dyy; if (trace * trace / det > MAX_CURVATURE_RATIO) continue; candidates.addElement(new float[] { fx, fy, fi }); //candidates.addElement( new float[]{ x, y, i } ); } } } } }