All Downloads are FREE. Search and download functionalities are using the official Maven repository.

boofcv.alg.feature.orientation.OrientationHistogramSift Maven / Gradle / Ivy

Go to download

BoofCV is an open source Java library for real-time computer vision and robotics applications.

There is a newer version: 1.1.7
Show newest version
/*
 * Copyright (c) 2011-2017, Peter Abeles. All Rights Reserved.
 *
 * This file is part of BoofCV (http://boofcv.org).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package boofcv.alg.feature.orientation;

import boofcv.alg.feature.detect.interest.FastHessianFeatureDetector;
import boofcv.core.image.FactoryGImageGray;
import boofcv.core.image.GImageGray;
import boofcv.misc.BoofMiscOps;
import boofcv.misc.CircularIndex;
import boofcv.numerics.InterpolateArray;
import boofcv.struct.ImageRectangle;
import boofcv.struct.image.ImageGray;
import georegression.metric.UtilAngle;
import org.ddogleg.struct.GrowQueue_F64;
import org.ddogleg.struct.GrowQueue_I32;

import java.util.Arrays;

/**
 * 

* Computes the orientation of a region around a point in scale-space as specified in the SIFT [1] paper. A histogram * of angles is computed using a weighted sum of image derivatives. The size of the region is specified by the * scale function parameter. Every pixel inside the sample region is read and contributes to the angle estimation. * If the image border is encountered the sample return is truncated. *

* *

* To get the orientation for the largest peak invoke {@link #getPeakOrientation()}. Other *

* * Differences from paper: *
    *
  • The angle in each bin is set to the atan2(y,x) of the weighted sum of image derivative
  • *
  • Interpolation is done using a 2nd degree polynomial instead of a parabola.
  • *
* *

* [1] Lowe, D. "Distinctive image features from scale-invariant keypoints". * International Journal of Computer Vision, 60, 2 (2004), pp.91--110. *

* * @author Peter Abeles */ public class OrientationHistogramSift { // How much does it inflate the scale by private double sigmaEnlarge; // Storage for orientation histogram. Each bin is for angles from i*histAngleBin to (i+1)*histAngleBin double histogramMag[]; // histograms containing the sum of each derivative double histogramX[]; double histogramY[]; // Number of radians each bin corresponds to. histAngleBin = 2*PI/histogram.length private double histAngleBin; // peaks in histogram private GrowQueue_I32 peaks = new GrowQueue_I32(10); // Found orientations of the feature private GrowQueue_F64 angles = new GrowQueue_F64(10); // Angle of the largest peak private double peakAngle; // local region from which the orientation is computed private ImageRectangle bound = new ImageRectangle(); // spacial image gradient of closest image in scale-space private GImageGray derivX,derivY; InterpolateArray approximateGauss; double approximateStep = 0.1; /** * Configures orientation estimation * * @param histogramSize Number of elements in the histogram. Standard is 36 * @param sigmaEnlarge How much the scale is enlarged by. Standard is 1.5 */ public OrientationHistogramSift(int histogramSize , double sigmaEnlarge , Class derivType ) { this.histogramMag = new double[ histogramSize ]; this.histogramX = new double[ histogramSize ]; this.histogramY = new double[ histogramSize ]; this.sigmaEnlarge = sigmaEnlarge; this.histAngleBin = 2.0*Math.PI/histogramSize; // compute an approximation of a Gaussian distribution as a function of the distance squared double samples[] = new double[ (int)(4*4/approximateStep) ]; for( int i = 0; i < samples.length; i++ ) { double dx2 = i*approximateStep; samples[i] = Math.exp(-0.5*dx2 ); } approximateGauss = new InterpolateArray(samples); this.derivX = FactoryGImageGray.create(derivType); this.derivY = FactoryGImageGray.create(derivType); } /** * Specify the input image */ public void setImageGradient(Deriv derivX, Deriv derivY ) { this.derivX.wrap(derivX); this.derivY.wrap(derivY); } /** * Estimates the orientation(s) of a region at the specified location and scale * * @param c_x Location x-axis * @param c_y Location y-axis * @param sigma blur standard deviations of detected feature. Also referred to as scale. */ public void process( double c_x , double c_y , double sigma ) { // convert to image coordinates int x = (int)(c_x + 0.5); int y = (int)(c_y + 0.5); // Estimate its orientation(s) computeHistogram(x, y, sigma ); // compute the descriptor findHistogramPeaks(); } /** * Constructs the histogram around the specified point. * * @param c_x Center x-axis * @param c_y Center y-axis * @param sigma Scale of feature, adjusted for local octave */ void computeHistogram(int c_x, int c_y, double sigma) { int r = (int)Math.ceil(sigma * sigmaEnlarge); // specify the area being sampled bound.x0 = c_x - r; bound.y0 = c_y - r; bound.x1 = c_x + r + 1; bound.y1 = c_y + r + 1; ImageGray rawDX = derivX.getImage(); ImageGray rawDY = derivY.getImage(); // make sure it is contained in the image bounds BoofMiscOps.boundRectangleInside(rawDX,bound); // clear the histogram Arrays.fill(histogramMag,0); Arrays.fill(histogramX,0); Arrays.fill(histogramY,0); // construct the histogram for( int y = bound.y0; y < bound.y1; y++ ) { // iterate through the raw array for speed int indexDX = rawDX.startIndex + y*rawDX.stride + bound.x0; int indexDY = rawDY.startIndex + y*rawDY.stride + bound.x0; for( int x = bound.x0; x < bound.x1; x++ ) { float dx = derivX.getF(indexDX++); float dy = derivY.getF(indexDY++); // edge intensity and angle double magnitude = Math.sqrt(dx*dx + dy*dy); double theta = UtilAngle.domain2PI(Math.atan2(dy,dx)); // weight from gaussian double weight = computeWeight( x-c_x, y-c_y , sigma ); // histogram index int h = (int)(theta / histAngleBin) % histogramMag.length; // update the histogram histogramMag[h] += magnitude*weight; histogramX[h] += dx*weight; histogramY[h] += dy*weight; } } } /** * Finds local peaks in histogram and selects orientations. Location of peaks is interpolated. */ void findHistogramPeaks() { // reset data structures peaks.reset(); angles.reset(); peakAngle = 0; // identify peaks and find the highest peak double largest = 0; int largestIndex = -1; double before = histogramMag[ histogramMag.length-2 ]; double current = histogramMag[ histogramMag.length-1 ]; for(int i = 0; i < histogramMag.length; i++ ) { double after = histogramMag[ i ]; if( current > before && current > after ) { int currentIndex = CircularIndex.addOffset(i,-1,histogramMag.length); peaks.push(currentIndex); if( current > largest ) { largest = current; largestIndex = currentIndex; } } before = current; current = after; } if( largestIndex < 0 ) return; // see if any of the other peaks are within 80% of the max peak double threshold = largest*0.8; for( int i = 0; i < peaks.size; i++ ) { int index = peaks.data[i]; current = histogramMag[index]; if( current >= threshold) { double angle = computeAngle(index); angles.push( angle ); if( index == largestIndex ) peakAngle = angle; } } } /** * Compute the angle. The angle for each neighbor bin is found using the weighted sum * of the derivative. Then the peak index is found by 2nd order polygon interpolation. These two bits of * information are combined and used to return the final angle output. * * @param index1 Histogram index of the peak * @return angle of the peak. -pi to pi */ double computeAngle( int index1 ) { int index0 = CircularIndex.addOffset(index1,-1, histogramMag.length); int index2 = CircularIndex.addOffset(index1, 1, histogramMag.length); // compute the peak location using a second order polygon double v0 = histogramMag[index0]; double v1 = histogramMag[index1]; double v2 = histogramMag[index2]; double offset = FastHessianFeatureDetector.polyPeak(v0,v1,v2); // interpolate using the index offset and angle of its neighbor return interpolateAngle(index0, index1, index2, offset); } /** * Given the interpolated index, compute the angle from the 3 indexes. The angle for each index * is computed from the weighted gradients. * @param offset Interpolated index offset relative to index0. range -1 to 1 * @return Interpolated angle. */ double interpolateAngle(int index0, int index1, int index2, double offset) { double angle1 = Math.atan2(histogramY[index1],histogramX[index1]); double deltaAngle; if( offset < 0 ) { double angle0 = Math.atan2(histogramY[index0],histogramX[index0]); deltaAngle = UtilAngle.dist(angle0,angle1); } else { double angle2 = Math.atan2(histogramY[index2], histogramX[index2]); deltaAngle = UtilAngle.dist(angle2,angle1); } return UtilAngle.bound(angle1 + deltaAngle*offset); } /** * Computes the weight based on a centered Gaussian shaped function. Interpolation is used to speed up the process */ double computeWeight( double deltaX , double deltaY , double sigma ) { // the exact equation // return Math.exp(-0.5 * ((deltaX * deltaX + deltaY * deltaY) / (sigma * sigma))); // approximation below. when validating this approach it produced results that were within // floating point tolerance of the exact solution, but much faster double d = ((deltaX * deltaX + deltaY * deltaY) / (sigma * sigma))/approximateStep; if( approximateGauss.interpolate(d) ) { return approximateGauss.value; } else return 0; } /** * A list of found orientations * * @return orientations */ public GrowQueue_F64 getOrientations() { return angles; } /** * Orientation of the largest peak */ public double getPeakOrientation() { return peakAngle; } }




© 2015 - 2025 Weber Informatics LLC | Privacy Policy