KeyPointDescriptor.cpp
Go to the documentation of this file.
00001 /*
00002 * This file is part of Parallel SURF, which implements the SURF algorithm
00003 * using multi-threading.
00004 *
00005 * Copyright (C) 2010 David Gossow
00006 *
00007 * It is based on the SURF implementation included in Pan-o-matic 0.9.4, 
00008 * written by Anael Orlinski.
00009 * 
00010 * Parallel SURF is free software; you can redistribute it and/or modify
00011 * it under the terms of the GNU General Public License as published by
00012 * the Free Software Foundation; either version 3 of the License, or
00013 * (at your option) any later version.
00014 * 
00015 * Parallel SURF is distributed in the hope that it will be useful,
00016 * but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 * GNU General Public License for more details.
00019 * 
00020 * You should have received a copy of the GNU General Public License
00021 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00022 */
00023 
00024 #include <algorithm>
00025 #include <vector>
00026 
00027 #include "KeyPointDescriptor.h"
00028 #include "KeyPoint.h"
00029 #include "MathStuff.h"
00030 #include "WaveFilter.h"
00031 
00032 using namespace parallelsurf;
00033 using namespace std;
00034 
00035 namespace parallelsurf
00036 {
00037   LUT<0, 83> Exp1(exp, 0.5, -0.08);
00038   LUT<0, 40> Exp2(exp, 0.5, -0.125);
00039 
00040   bool operator < ( const response a, const response b )
00041   {
00042     return a.orientation < b.orientation;
00043   }
00044 }
00045 
00046 KeyPointDescriptor::KeyPointDescriptor(Image& iImage, boost::threadpool::pool &iThreadPool, bool iExtended) : 
00047   _image(iImage), _extended(iExtended), _threadPool(iThreadPool)
00048 {
00049   // init some parameters
00050   _subRegions = 4;
00051 
00052   _vecLen = 4;
00053   if (_extended)
00054     _vecLen = 8;
00055   
00056   _magFactor = 12.0 / _subRegions;
00057 }
00058 
00059 void KeyPointDescriptor::makeDescriptor(KeyPoint& ioKeyPoint) const
00060 {
00061 //   std::cout << "x";
00062   // create a descriptor context
00063   KeyPointDescriptorContext aCtx(_subRegions, _vecLen, ioKeyPoint._ori);
00064 
00065   // create the storage in the keypoint
00066   ioKeyPoint._vec.resize(getDescriptorLength());
00067 
00068   // assign the orientation
00069   //assignOrientation(ioKeyPoint);
00070 
00071   // create a vector
00072   createDescriptor(aCtx, ioKeyPoint);
00073 
00074   // transform back to vector
00075 
00076   // fill the vector with the values of the square...
00077   // remember the shift of 1 to drop outborders.
00078   size_t aV = 0;
00079   for(int aYIt = 1; aYIt < _subRegions+1; ++aYIt)
00080   {
00081     for(int aXIt = 1; aXIt < _subRegions+1; ++aXIt)
00082     {
00083       for(int aVIt = 0; aVIt < _vecLen; ++aVIt)
00084       {
00085         double a = aCtx._cmp[aYIt][aXIt][aVIt];
00086         ioKeyPoint._vec[aV] = a;
00087         aV++;
00088       }
00089     }
00090   }
00091 
00092   // normalize
00093   Math::Normalize(ioKeyPoint._vec);
00094 }
00095 
00096 
00097 void KeyPointDescriptor::assignOrientation(KeyPoint& ioKeyPoint) const
00098 {
00099   unsigned int aRX = Math::Round ( ioKeyPoint._x );
00100   unsigned int aRY = Math::Round ( ioKeyPoint._y );
00101   int aStep = ( int ) ( ioKeyPoint._scale + 0.8 );
00102   
00103   WaveFilter aWaveFilter ( 2.0 * ioKeyPoint._scale + 1.6, _image );
00104   
00105   vector< response > aRespVector;
00106   aRespVector.reserve ( 253 );
00107   
00108   // compute haar wavelet responses in a circular neighborhood of 6s
00109   for ( int aYIt = -9; aYIt <= 9; aYIt++ )
00110   {
00111     int aSY = aRY + aYIt * aStep;
00112     for ( int aXIt = -9; aXIt <= 9; aXIt++ )
00113     {
00114       int aSX = aRX + aXIt * aStep;
00115       
00116       // keep points in a circular region of diameter 6s
00117       unsigned int aSqDist = aXIt * aXIt + aYIt * aYIt;
00118       if ( aSqDist <= 81 && aWaveFilter.checkBounds ( aSX, aSY ) )
00119       {
00120         double aWavX = aWaveFilter.getWx ( aSX, aSY );
00121         double aWavY = aWaveFilter.getWy ( aSX, aSY );
00122         double aWavResp = sqrt ( aWavX * aWavX + aWavY * aWavY );
00123         if ( aWavResp > 0 )
00124         {
00125           double aAngle = atan2 ( aWavY, aWavX );
00126           response resp;
00127           resp.orientation =  aAngle;
00128           resp.magnitude = aWavResp * Exp1 ( aSqDist );
00129           aRespVector.push_back ( resp );
00130         }
00131       }
00132     }
00133   }
00134   
00135   //assert ( aRespVector.size() <= 253 );
00136   
00137   //no wavelet responses != 0, can't calculate orientation
00138   if ( aRespVector.size() == 0 )
00139   {
00140     ioKeyPoint._ori = 0;
00141     return;
00142   }
00143   
00144   //sort responses by orientation
00145   sort ( aRespVector.begin(), aRespVector.end() );
00146   
00147   //estimate orientation using a sliding window of PI/3
00148   response aMax = aRespVector[0];
00149   double aMagnitudeSum = aRespVector[0].magnitude;
00150   double aOrientationSum = aRespVector[0].orientation * aMagnitudeSum;
00151   
00152   size_t aWindowBegin = 0;
00153   size_t aWindowEnd = 0;
00154   
00155   float aOriAdd = 0;
00156   
00157   while ( aWindowBegin < aRespVector.size() )
00158   {
00159     float aWindowSize = (aRespVector[aWindowEnd].orientation + aOriAdd) - aRespVector[aWindowBegin].orientation;
00160     if ( aWindowSize < PI / 3 )
00161     {
00162       //found new max.
00163       if ( aMagnitudeSum > aMax.magnitude )
00164       {
00165         aMax.magnitude = aMagnitudeSum;
00166         aMax.orientation = aOrientationSum;
00167       }
00168       aWindowEnd++;
00169       if ( aWindowEnd >= aRespVector.size() )
00170       {
00171         aWindowEnd = 0;
00172         aOriAdd += 2*PI;
00173       }
00174       aMagnitudeSum += aRespVector[aWindowEnd].magnitude;
00175       aOrientationSum += aRespVector[aWindowEnd].magnitude * (aRespVector[aWindowEnd].orientation + aOriAdd);
00176     }
00177     else
00178     {
00179       aMagnitudeSum -= aRespVector[aWindowBegin].magnitude;
00180       aOrientationSum -= aRespVector[aWindowBegin].magnitude * (aRespVector[aWindowEnd].orientation + aOriAdd);
00181       aWindowBegin++;
00182     }
00183   }
00184   
00185   ioKeyPoint._ori = aMax.orientation / aMax.magnitude;
00186 }
00187 
00188 
00189 void KeyPointDescriptor::createDescriptor(KeyPointDescriptorContext& iCtx, KeyPoint& ioKeyPoint) const
00190 {
00191   // create the vector of features by analyzing a square patch around the point.
00192   // for this the current patch (x,y) will be translated in rotated coordinates (u,v)
00193 
00194   double aX = ioKeyPoint._x;
00195   double aY = ioKeyPoint._y;
00196   double aS = ioKeyPoint._scale * 1.65; // multiply by this nice constant
00197 
00198   // make integer values from double ones
00199   int aIntX = Math::Round(aX);
00200   int aIntY = Math::Round(aY);
00201   int aIntS = Math::Round(aS / 2.0);  
00202   if (aIntS < 1) aIntS = 1;
00203 
00204   // calc subpixel shift
00205   double aSubX = aX - aIntX;
00206   double aSubY = aY - aIntY;
00207 
00208   // calc subpixel shift in rotated coordinates
00209   double aSubV = iCtx._cos * aSubY + iCtx._sin * aSubX;
00210   double aSubU = - iCtx._sin * aSubY + iCtx._cos * aSubX;
00211 
00212   // calc step of sampling
00213   double aStepSample = aS * _magFactor;
00214 
00215   // make a bounding box around the rotated patch square. 
00216   double aRadius = (1.414 * aStepSample) * (_subRegions + 1) / 2.0;
00217   int aIntRadius = Math::Round(aRadius / aIntS);
00218   if (aS < 1) aS = 1;
00219 
00220   // go through all the pixels in the bounding box
00221   for (int aYIt = -aIntRadius; aYIt <= aIntRadius; ++aYIt)
00222   {
00223     for (int aXIt = -aIntRadius; aXIt <= aIntRadius; ++aXIt)
00224     {
00225       // calculate U,V rotated values from X,Y taking in account subpixel correction
00226       // divide by step sample to put in index units
00227       double aU = ((( - iCtx._sin * aYIt + iCtx._cos * aXIt) * aIntS) - aSubU) / aStepSample;
00228       double aV = (((iCtx._cos * aYIt + iCtx._sin * aXIt) * aIntS) - aSubV) / aStepSample;
00229 
00230       // compute location of sample in terms of real array coordinates 
00231       double aUIdx = _subRegions / 2.0 - 0.5 + aU;
00232       double aVIdx = _subRegions / 2.0 - 0.5 + aV;
00233 
00234       // test if some bins will be filled.
00235       if (aUIdx >= -1.0 && aUIdx < _subRegions &&
00236         aVIdx >= -1.0 && aVIdx < _subRegions )
00237       {
00238         int aXSample = aIntS * aXIt + aIntX;
00239         int aYSample = aIntS * aYIt + aIntY;
00240         int aStep = (int)aS;
00241 
00242         WaveFilter aWaveFilter(aStep, _image);
00243         if (!aWaveFilter.checkBounds(aXSample, aYSample))
00244           continue;
00245 
00246         double aExp = Exp2((int)(aV * aV + aU * aU)); 
00247 
00248         double aWavX = aWaveFilter.getWx(aXSample, aYSample) * aExp;
00249         double aWavY = aWaveFilter.getWy(aXSample, aYSample) * aExp;
00250 
00251         double aWavXR = (iCtx._cos * aWavX + iCtx._sin * aWavY);
00252         double aWavYR = (iCtx._sin * aWavX - iCtx._cos * aWavY);
00253 
00254         // due to the rotation, the patch has to be dispatched in 2 bins in each direction
00255         // get the bins and weight for each of them
00256         // shift by 1 to avoid checking bounds
00257         const int aBin1U = (int)(aUIdx + 1.0);
00258         const int aBin2U = aBin1U + 1;
00259         const int aBin1V = (int)(aVIdx + 1.0);
00260         const int aBin2V = aBin1V + 1;
00261 
00262         const double aWeightBin1U = aBin1U - aUIdx;
00263         const double aWeightBin2U = 1 - aWeightBin1U;
00264 
00265         const double aWeightBin1V = aBin1V - aVIdx;
00266         const double aWeightBin2V = 1 - aWeightBin1V;
00267 
00268         if(_extended) 
00269         {
00270           int aBin = (aWavYR <= 0) ? 0 : 1;
00271           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aWavXR;
00272           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aWavXR;
00273           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aWavXR;
00274           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aWavXR;
00275 
00276           aBin += 2;
00277           double aVal = fabs(aWavXR);
00278           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aVal;
00279           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aVal;
00280           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aVal;
00281           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aVal;
00282 
00283           aBin = (aWavXR <= 0) ? 4 : 5;
00284           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aWavYR;
00285           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aWavYR;
00286           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aWavYR;
00287           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aWavYR;
00288 
00289           aBin += 2;
00290           aVal = fabs(aWavYR);
00291           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aVal;
00292           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aVal;
00293           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aVal;
00294           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aVal;
00295 
00296         } 
00297         else 
00298         {
00299           int aBin = (aWavXR <= 0) ? 0 : 1;
00300           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aWavXR;
00301           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aWavXR;
00302           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aWavXR;
00303           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aWavXR;
00304 
00305           aBin = (aWavYR <= 0) ? 2 : 3;
00306           iCtx._cmp[aBin1V][aBin1U][aBin] += aWeightBin1V * aWeightBin1U * aWavYR;
00307           iCtx._cmp[aBin2V][aBin1U][aBin] += aWeightBin2V * aWeightBin1U * aWavYR;
00308           iCtx._cmp[aBin1V][aBin2U][aBin] += aWeightBin1V * aWeightBin2U * aWavYR;
00309           iCtx._cmp[aBin2V][aBin2U][aBin] += aWeightBin2V * aWeightBin2U * aWavYR;
00310 
00311         }
00312       }   
00313     }
00314   }
00315 }
00316 
00317 
00318 KeyPointDescriptorContext::KeyPointDescriptorContext( int iSubRegions, 
00319                             int iVecLen, 
00320                             double iOrientation) :
00321   _subRegions(iSubRegions), _sin(sin(iOrientation)), _cos(cos(iOrientation))
00322 {
00323   // allocate and initialize the components of the vector
00324   // the idea is to allocate 2 more in each direction and
00325   // shift access by 1 to discard out of bounds.
00326 
00327   int aExtSub = _subRegions + 2;
00328   _cmp = new double **[aExtSub];
00329   for (int aYIt = 0; aYIt < aExtSub; ++aYIt)
00330   {
00331     _cmp[aYIt] = new double *[aExtSub];
00332     for (int aXIt = 0; aXIt < aExtSub; ++aXIt)
00333     {
00334       _cmp[aYIt][aXIt] = new double[iVecLen];
00335       for (int aVIt = 0; aVIt < iVecLen; ++aVIt)
00336         _cmp[aYIt][aXIt][aVIt] = 0;
00337     }
00338   }
00339 }
00340 
00341 
00342 KeyPointDescriptorContext::~KeyPointDescriptorContext()
00343 {
00344   // destroy the components of the vector
00345   int aExtSub = _subRegions + 2;
00346   for (int aYIt = 0; aYIt < aExtSub; ++aYIt)
00347   {
00348     for (int aXIt = 0; aXIt < aExtSub; ++aXIt)
00349     {
00350       delete[] _cmp[aYIt][aXIt];
00351     }
00352     delete[] _cmp[aYIt];
00353   }
00354 }
00355 
00356 
00357 int KeyPointDescriptor::getDescriptorLength() const
00358 {
00359   return _vecLen * _subRegions * _subRegions;
00360 }


or_libs
Author(s): raphael
autogenerated on Mon Oct 6 2014 02:53:18