001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.image.feature.local.detector.pyramid;
031
032import org.openimaj.citation.annotation.Reference;
033import org.openimaj.citation.annotation.ReferenceType;
034import org.openimaj.citation.annotation.References;
035import org.openimaj.image.FImage;
036import org.openimaj.image.analysis.pyramid.gaussian.GaussianOctave;
037
038import Jama.Matrix;
039
040/**
041 * Implementation of the method described in
042 * "Invariant Features from Interest Point Groups" by Matthew Brown and David
043 * Lowe (http://www.cs.ubc.ca/~lowe/papers/brown02.pdf) for improving the
044 * localisation of interest points detected in a difference-of-Gaussian by
045 * fitting a 3D quadratic to the scale-space Laplacian (approximated by the
046 * difference-of-Gaussian pyramid).
047 * 
048 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
049 * 
050 */
051@References(references = {
052                @Reference(
053                                type = ReferenceType.Inproceedings,
054                                author = { "Matthew Brown", "D Lowe" },
055                                title = "Invariant Features from Interest Point Groups",
056                                year = "2002",
057                                booktitle = "BMVC 2002: 13th British Machine Vision Conference",
058                                pages = { "253", "", "262" },
059                                month = "September"
060                ),
061                @Reference(
062                                type = ReferenceType.Article,
063                                author = { "David Lowe" },
064                                title = "Distinctive image features from scale-invariant keypoints",
065                                year = "2004",
066                                journal = "IJCV",
067                                pages = { "91", "110" },
068                                month = "January",
069                                number = "2",
070                                volume = "60")
071})
072public class InterpolatingOctaveExtremaFinder extends BasicOctaveExtremaFinder {
073        /**
074         * Default number of interpolation iterations
075         */
076        public static final int DEFAULT_INTERPOLATION_ITERATIONS = 5;
077
078        protected int numInterpolationIterations;
079
080        private boolean[][] map;
081        private int currentIteration;
082
083        /**
084         * Default constructor using {@link #DEFAULT_MAGNITUDE_THRESHOLD} for the
085         * magnitude threshold, {@link #DEFAULT_EIGENVALUE_RATIO} for the Eigenvalue
086         * ratio threshold and {@link #DEFAULT_INTERPOLATION_ITERATIONS} for the
087         * number of iterations.
088         */
089        public InterpolatingOctaveExtremaFinder() {
090                this(DEFAULT_MAGNITUDE_THRESHOLD, DEFAULT_EIGENVALUE_RATIO, DEFAULT_INTERPOLATION_ITERATIONS);
091        }
092
093        /**
094         * Construct with the given magnitude threshold,
095         * {@link #DEFAULT_EIGENVALUE_RATIO} for the Eigenvalue ratio threshold and
096         * {@link #DEFAULT_INTERPOLATION_ITERATIONS} for the number of iterations.
097         * 
098         * @param magnitudeThreshold
099         *            the magnitude threshold
100         */
101        public InterpolatingOctaveExtremaFinder(float magnitudeThreshold) {
102                this(magnitudeThreshold, DEFAULT_EIGENVALUE_RATIO, DEFAULT_INTERPOLATION_ITERATIONS);
103        }
104
105        /**
106         * Construct with the given magnitude and Eigenvalue thresholds andnumber of
107         * iterations.
108         * 
109         * @param magnitudeThreshold
110         *            the magnitude threshold
111         * @param eigenvalueRatio
112         *            the Eigenvalue threshold
113         * @param numInterpolationIterations
114         *            the number of interpolation iterations
115         */
116        public InterpolatingOctaveExtremaFinder(float magnitudeThreshold, float eigenvalueRatio,
117                        int numInterpolationIterations)
118        {
119                super(magnitudeThreshold, eigenvalueRatio);
120                this.numInterpolationIterations = numInterpolationIterations;
121        }
122
123        @Override
124        public void process(GaussianOctave<FImage> octave) {
125                map = new boolean[octave.images[0].height][octave.images[0].width];
126
127                super.process(octave);
128        }
129
130        @Override
131        protected void processExtrema(FImage[] dogs, int s, int x, int y, float octSize) {
132                currentIteration = 0;
133                processExtremaInternal(dogs, s, x, y, octSize);
134        }
135
136        protected void processExtremaInternal(FImage[] dogs, int s, int x, int y, float octSize) {
137                int newy = y, newx = x;
138
139                // fit 3d quadratic
140                final FitResult fit = fitQuadratic3D(dogs, s, x, y);
141
142                if (fit.offset.get(1, 0) > 0.5 && y < dogs[0].height - octave.options.getBorderPixels())
143                        newy++;
144                if (fit.offset.get(1, 0) < -0.5 && y > octave.options.getBorderPixels())
145                        newy--;
146                if (fit.offset.get(2, 0) > 0.5 && x < dogs[0].width - octave.options.getBorderPixels())
147                        newx++;
148                if (fit.offset.get(2, 0) < -0.5 && x > octave.options.getBorderPixels())
149                        newx--;
150                if (currentIteration < numInterpolationIterations && (newy != y || newx != x)) {
151                        currentIteration++;
152                        processExtremaInternal(dogs, s, newx, newy, octSize);
153                        return;
154                }
155
156                // if the offset is too great, or the peak doesn't have enough contrast
157                // w.r.t. the threshold, then stop
158                if (Math.abs(fit.offset.get(0, 0)) > 1.5 || Math.abs(fit.offset.get(1, 0)) > 1.5
159                                || Math.abs(fit.offset.get(2, 0)) > 1.5 || Math.abs(fit.peakval) < normMagnitudeScales)
160                        return;
161
162                // avoid duplicates by checking a map of locations
163                if (map[y][x] == true)
164                        return;
165                map[y][x] = true;
166
167                // calculate the new scale
168                final float octaveScale = octave.options.getInitialSigma()
169                                * (float) Math.pow(2.0, (s + fit.offset.get(0, 0)) / octave.options.getScales());
170
171                // fire the listener
172                if (listener != null)
173                        listener.foundInterestPoint(this, x + (float) fit.offset.get(2, 0), y + (float) fit.offset.get(1, 0),
174                                        octaveScale);
175        }
176
177        class FitResult {
178                Matrix offset; // the offset of 0 or the peak
179                float peakval; // the interpolated value at the peak
180        }
181
182        // fit a quadratic around the point
183        FitResult fitQuadratic3D(FImage[] dogs, int s, int x, int y) {
184                final float[][] dog0 = dogs[s - 1].pixels;
185                final float[][] dog1 = dogs[s].pixels;
186                final float[][] dog2 = dogs[s + 1].pixels;
187
188                final Matrix H = new Matrix(3, 3); // hessian matrix; derivatives
189                                                                                        // estimated by differences across
190                                                                                        // scales
191                H.set(0, 0, dog0[y][x] - 2.0 * dog1[y][x] + dog2[y][x]);
192                H.set(0, 1, ((dog2[y + 1][x] - dog2[y - 1][x]) - (dog0[y + 1][x] - dog0[y - 1][x])) / 4.0);
193                H.set(0, 2, ((dog2[y][x + 1] - dog2[y][x - 1]) - (dog0[y][x + 1] - dog0[y][x - 1])) / 4.0);
194
195                H.set(1, 0, H.get(0, 1));
196                H.set(1, 1, dog1[y - 1][x] - 2.0 * dog1[y][x] + dog1[y + 1][x]);
197                H.set(1, 2, ((dog1[y + 1][x + 1] - dog1[y + 1][x - 1]) - (dog1[y - 1][x + 1] - dog1[y - 1][x - 1])) / 4.0);
198
199                H.set(2, 0, H.get(0, 2));
200                H.set(2, 1, H.get(1, 2));
201                H.set(2, 2, dog1[y][x - 1] - 2.0 * dog1[y][x] + dog1[y][x + 1]);
202
203                // gradient vector
204                final Matrix gM = new Matrix(
205                                new double[][] {
206                                                { (dog2[y][x] - dog0[y][x]) / 2.0f }, { (dog1[y + 1][x] - dog1[y - 1][x]) / 2.0f }, { (dog1[y][x + 1] - dog1[y][x - 1]) / 2.0f }
207                                });
208
209                final Matrix offsetM = H.solve(gM.times(-1));
210
211                final FitResult result = new FitResult();
212                result.offset = offsetM;
213
214                final float dp = (float) (offsetM.get(0, 0) * gM.get(0, 0) + offsetM.get(1, 0) * gM.get(1, 0) + offsetM.get(2, 0)
215                                * gM.get(2, 0));
216                result.peakval = (dog1[y][x] + 0.5f * dp);
217
218                return result;
219        }
220}