CircularKeyPointDescriptor.cpp

Go to the documentation of this file.
00001 /* *-* tab-width: 4; c-basic-offset: 4; -*- */
00002 /*
00003 * Copyright (C) 2007-2008 Anael Orlinski
00004 *
00005 * This file is part of Panomatic.
00006 *
00007 * Panomatic is free software; you can redistribute it and/or modify
00008 * it under the terms of the GNU General Public License as published by
00009 * the Free Software Foundation; either version 2 of the License, or
00010 * (at your option) any later version.
00011 *
00012 * Panomatic is distributed in the hope that it will be useful,
00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 * GNU General Public License for more details.
00016 *
00017 * You should have received a copy of the GNU General Public License
00018 * along with Panomatic; if not, write to the Free Software
00019 * <http://www.gnu.org/licenses/>.
00020 */
00021 
00022 #include <iostream>
00023 
00024 #ifdef _MSC_VER
00025 #define _USE_MATH_DEFINES
00026 #endif
00027 #include <math.h>
00028 #include <vector>
00029 #include <map>
00030 #include <fstream>
00031 #include <cassert>
00032 
00033 #include <string.h>
00034 
00035 #include "KeyPoint.h"
00036 #include "CircularKeyPointDescriptor.h"
00037 #include "MathStuff.h"
00038 #include "WaveFilter.h"
00039 
00040 //#define DEBUG_DESC
00041 //#define DEBUG_ROT
00042 
00043 using namespace lfeat;
00044 using namespace std;
00045 
00046 LUT<0, 83> Exp1_2(exp, 0.5, -0.08);
00047 
00048 CircularKeyPointDescriptor::CircularKeyPointDescriptor(Image& iImage,
00049         std::vector<int> rings, std::vector<double>ring_radius,
00050         std::vector<double>ring_gradient_width,
00051         int ori_nbins, double ori_sample_scale, int ori_gridsize) :
00052     _image(iImage), _ori_nbins(ori_nbins), _ori_sample_scale(ori_sample_scale),
00053     _ori_gridsize(ori_gridsize)
00054 {
00055     // default parameters
00056     if (rings.size() == 0)
00057     {
00058         rings.push_back(1);
00059         ring_radius.push_back(0);
00060         ring_gradient_width.push_back(2);
00061         rings.push_back(6);
00062         ring_radius.push_back(4.2);
00063         ring_gradient_width.push_back(1.9);
00064         /*
00065         rings.push_back(6);
00066         ring_radius.push_back(5);
00067         ring_gradient_width.push_back(3);
00068         */
00069         rings.push_back(6);
00070         ring_radius.push_back(8);
00071         ring_gradient_width.push_back(3.2);
00072     }
00073 
00074     assert(rings.size() == ring_radius.size());
00075     assert(rings.size() == ring_gradient_width.size());
00076 
00077     /*
00078     // radius of the rings (ring 0 is the center point)
00079     double ring_radius[3];
00080     ring_radius[0] =  0;
00081     ring_radius[1] =  2;
00082     ring_radius[2] =  5;
00083 
00084     // width of the boxfilter used during the computation (ring 0 is the center point)
00085     double ring_gradwidth[3];
00086     ring_gradwidth[0] = 0.5;
00087     ring_gradwidth[1] = 1.5;
00088     ring_gradwidth[2] = 2.5;
00089 
00090     // number of entries for each ring
00091     int nrings[3];
00092     nrings[0] = 1;
00093     nrings[1] = 6;
00094     nrings[2] = 6;
00095     */
00096     // compute number of sampling points
00097     _subRegions = 0;
00098     for (unsigned int i=0; i < rings.size(); i++)
00099     {
00100         _subRegions += rings[i];
00101     }
00102 
00103     // create a list of sample parameters
00104     _samples = new SampleSpec[_subRegions];
00105 
00106     // precompute positions of the sampling points
00107     int j=0;
00108     for (unsigned int i=0; i < rings.size(); i++)
00109     {
00110         // alternate the phi offset of the rings,
00111         // so that the circles don't overlap too much with the
00112         // next ring
00113         double phioffset = i % 2== 0 ? 0 : M_PI/rings[i];
00114         for (int ri=0; ri < rings[i]; ri++)
00115         {
00116             double phi = ri*2*M_PI /rings[i] + phioffset;
00117             _samples[j].x = ring_radius[i]*cos(phi);
00118             _samples[j].y = ring_radius[i]*sin(phi);
00119             _samples[j].size = ring_gradient_width[i];
00120             j++;
00121         }
00122     }
00123 
00124     // compute 4 gradient entries (pos(dx), pos(-dx), pos(dy), pos(-dy))
00125     // maybe use something else here, for example, just the gradient entries...
00126     // also use LBP as a 5th feature
00127     _vecLen = 3;
00128     _descrLen = _vecLen * _subRegions - 1;
00129 
00130     _ori_hist = new double [_ori_nbins + 2];
00131 
00132 }
00133 
00134 CircularKeyPointDescriptor::~CircularKeyPointDescriptor()
00135 {
00136     delete[] _ori_hist;
00137     delete[] _samples;
00138 }
00139 
00140 void CircularKeyPointDescriptor::makeDescriptor(KeyPoint& ioKeyPoint) const
00141 {
00142     // create a descriptor context
00143     //KeyPointDescriptorContext aCtx(_subRegions, _vecLen, ioKeyPoint._ori);
00144 
00145     // create the storage in the keypoint
00146     if (!ioKeyPoint._vec)
00147     {
00148         ioKeyPoint.allocVector(getDescriptorLength());
00149     }
00150 
00151     // create a vector
00152     createDescriptor(ioKeyPoint);
00153 
00154     // normalize
00155     Math::Normalize(ioKeyPoint._vec, getDescriptorLength());
00156 }
00157 
00158 int CircularKeyPointDescriptor::assignOrientation(KeyPoint& ioKeyPoint, double angles[4]) const
00159 {
00160     double* hist = _ori_hist+1;
00161     unsigned int aRX = Math::Round(ioKeyPoint._x);
00162     unsigned int aRY = Math::Round(ioKeyPoint._y);
00163     int aStep = (int)(ioKeyPoint._scale + 0.8);
00164 
00165     WaveFilter aWaveFilter(_ori_sample_scale * ioKeyPoint._scale + 1.5, _image);
00166 
00167 #ifdef DEBUG_ROT_2
00168     std::cerr << "ori_scale = " << 2.5 * ioKeyPoint._scale + 1.5 << std::endl;
00169     std::cerr << "ori= [ ";
00170 #endif
00171 
00172 #ifdef _MSC_VER
00173 #pragma message("use LUT after parameter tuning")
00174 #else
00175 #warning use LUT after parameter tuning
00176 #endif
00177     double coeffadd = 0.5;
00178     double coeffmul = (0.5 + 6 ) / - (_ori_nbins*_ori_nbins);
00179 
00180     memset(_ori_hist, 0, sizeof(double)*(_ori_nbins + 2));
00181     // compute haar wavelet responses in a circular neighborhood of _ori_gridsize s
00182     for (int aYIt = -_ori_gridsize; aYIt <= _ori_gridsize; aYIt++)
00183     {
00184         int aSY = aRY + aYIt * aStep;
00185         for (int aXIt = -_ori_gridsize; aXIt <= _ori_gridsize; aXIt++)
00186         {
00187             int aSX = aRX + aXIt * aStep;
00188             // keep points in a circular region of diameter 6s
00189             const int aSqDist = aXIt * aXIt + aYIt * aYIt;
00190             if (aSqDist <= _ori_nbins*_ori_nbins && aWaveFilter.checkBounds(aSX, aSY))
00191             {
00192                 double aWavX = aWaveFilter.getWx(aSX, aSY);
00193                 // y axis for derivate seems is from bottom to top (typical math coordinate system),
00194                 // not from top to bottom (as in the typical image coordinate system).
00195                 double aWavY = - aWaveFilter.getWy(aSX, aSY);
00196                 double aWavResp = sqrt(aWavX * aWavX + aWavY * aWavY);
00197                 if (aWavResp > 0)
00198                 {
00199                     // add PI ->  0 .. 2*PI interval
00200                     double angle = atan2(aWavY, aWavX)+PI;
00201                     int bin =  angle/(2*PI) * _ori_nbins;
00202                     // deal with possible rounding problems.
00203                     bin = (bin+_ori_nbins)%_ori_nbins;
00204                     // center of bin 0 equals -PI + 16┬░deg, etc.
00205                     double weight = exp(coeffmul * (aSqDist+coeffadd));
00206                     hist[bin] += aWavResp * weight;
00207                     //hist[bin] += aWavResp * Exp1_2(aSqDist);
00208 #ifdef DEBUG_ROT_2
00209                     std::cerr << "[ " << aSX << ", " << aSY << ", "
00210                               << aWavX << ", " << aWavY << ", "
00211                               << aWavResp << ", " << angle << ", " << bin << ", "
00212                               << aSqDist << ", " << aWavResp* Exp1_2(aSqDist) << "], " << std::endl;
00213 #endif
00214                 }
00215             }
00216         }
00217     }
00218 #ifdef DEBUG_ROT_2
00219     std::cerr << "]" << endl;
00220 #endif
00221 
00222 #if 0
00223     // smoothing doesn't seem to work very well with the default grid + scale values for orientation histogram building
00224     // smooth histogram
00225     for (int it=0; it < 1; it++)
00226     {
00227         double prev = hist[_ori_nbins-1];
00228         double first = hist[0];
00229         for (int i=0; i < _ori_nbins-2; i++)
00230         {
00231             double hs = (prev + hist[i] + hist[i+1]) / 3.0;
00232             prev = hist[i];
00233             hist[i] = hs;
00234         }
00235         hist[_ori_nbins-1] = (prev + hist[_ori_nbins-1] + first) / 3.0;
00236     }
00237 #endif
00238 
00239     // avoid boundary problems, wrap around histogram
00240     hist[-1] = hist[_ori_nbins-1];
00241     hist[_ori_nbins] = hist[0];
00242 
00243     // find bin with the maximum response.
00244     double aMax = hist[0];
00245     int iMax = 0;
00246 #ifdef DEBUG_ROT
00247     std::cerr << "rot_hist: [ " << aMax;
00248 #endif
00249     for (int i=1; i < _ori_nbins; i++)
00250     {
00251 #ifdef DEBUG_ROT
00252         std::cerr << ", " << hist[i];
00253 #endif
00254         if (hist[i] > aMax)
00255         {
00256             aMax = hist[i];
00257             iMax = i;
00258         }
00259     }
00260 #ifdef DEBUG_ROT
00261     std::cerr << " ] " << std::endl;
00262 #endif
00263 
00264     double prev = hist[iMax-1];
00265     double curr = hist[iMax];
00266     double next = hist[iMax+1];
00267     double dsub = -0.5*(next-prev)/(prev+next-2*curr);
00268     ioKeyPoint._ori = (iMax+0.5+dsub) / _ori_nbins * 2*PI - PI;
00269 
00270     // add more keypoints that are within 0.8 of the maximum strength.
00271     aMax *= 0.8;
00272     int nNewOri = 0;
00273     for (int i=0; i < _ori_nbins; i++)
00274     {
00275         double prev = hist[i-1];
00276         double curr = hist[i];
00277         double next = hist[i+1];
00278         if (curr > aMax && prev < curr && next < curr && i != iMax)
00279         {
00280             dsub = -0.5*(next-prev)/(prev+next-2*curr);
00281             angles[nNewOri] = (i+0.5+dsub) / _ori_nbins * 2*PI - PI;
00282             nNewOri++;
00283             if (nNewOri == 4)
00284             {
00285                 break;
00286             }
00287         }
00288     }
00289     return nNewOri;
00290 }
00291 
00292 // gradient and intensity difference
00293 void CircularKeyPointDescriptor::createDescriptor(KeyPoint& ioKeyPoint) const
00294 {
00295 #ifdef DEBUG_DESC
00296     std::ofstream dlog("descriptor_details.txt", std::ios_base::app);
00297 #endif
00298 
00299     // create the vector of features by analyzing a square patch around the point.
00300     // for this the current patch (x,y) will be translated in rotated coordinates (u,v)
00301 
00302     double aX = ioKeyPoint._x;
00303     double aY = ioKeyPoint._y;
00304     int aS = (int)ioKeyPoint._scale;
00305 
00306     // get the sin/cos of the orientation
00307     double ori_sin = sin(ioKeyPoint._ori);
00308     double ori_cos = cos(ioKeyPoint._ori);
00309 
00310     if (aS < 1)
00311     {
00312         aS = 1;
00313     }
00314 
00315     // compute the gradients in x and y for all regions of interest
00316 
00317     // we override the wave filter size later
00318     WaveFilter aWaveFilter(10, _image);
00319 
00320     // compute features at each position and store in feature vector
00321     int j=0;
00322     double middleMean = 0;
00323 
00324     for (int i=0; i < _subRegions; i++)
00325     {
00326         // scale radius with aS.
00327         double xS = _samples[i].x * aS;
00328         double yS = _samples[i].y * aS;
00329         // rotate sample point with the orientation
00330         double aXSample = aX + xS * ori_cos - yS * ori_sin;
00331         double aYSample = aY + xS * ori_sin + yS * ori_cos;
00332         // make integer values from double ones
00333         int aIntXSample = Math::Round(aXSample);
00334         int aIntYSample = Math::Round(aYSample);
00335         int aIntSampleSize = Math::Round(_samples[i].size* aS);
00336         int sampleArea = aIntSampleSize * aIntSampleSize;
00337 
00338         if (!aWaveFilter.checkBounds(aIntXSample, aIntYSample, aIntSampleSize))
00339         {
00340             ioKeyPoint._vec[j++] = 0;
00341             ioKeyPoint._vec[j++] = 0;
00342             //ioKeyPoint._vec[j++] = 0;
00343             //ioKeyPoint._vec[j++] = 0;
00344             if (i > 0)
00345             {
00346                 ioKeyPoint._vec[j++] = 0;
00347             }
00348 #ifdef DEBUG_DESC
00349             dlog << xS << " " << yS << " "
00350                  << aIntXSample << " " << aIntYSample << " " << aIntSampleSize << " "
00351                  << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << std::endl;
00352 #endif
00353             continue;
00354         }
00355 
00356         double aWavX = aWaveFilter.getWx(aIntXSample, aIntYSample, aIntSampleSize) / sampleArea;
00357         double aWavY = - aWaveFilter.getWy(aIntXSample, aIntYSample, aIntSampleSize) / sampleArea;
00358 
00359         double meanGray = aWaveFilter.getSum(aIntXSample, aIntYSample, aIntSampleSize) / sampleArea;
00360 
00361         if (i == 0)
00362         {
00363             middleMean = meanGray;
00364         }
00365 
00366         // rotate extracted gradients
00367 
00368         //  need to rotate in the other direction?
00369 
00370         double aWavXR = aWavX * ori_cos + aWavY * ori_sin;
00371         double aWavYR = -aWavX * ori_sin + aWavY * ori_cos;
00372         /*
00373         double aWavXR = aWavX * ori_cos - aWavY * ori_sin;
00374         double aWavYR = aWavX * ori_sin + aWavY * ori_cos;
00375         */
00376 
00377 #ifdef DEBUG_DESC
00378         dlog << xS << " " << yS << " "
00379              << aIntXSample << " " << aIntYSample << " " << aIntSampleSize << " "
00380              << aWavX << " " << aWavY << " " << meanGray << " " << aWavXR << " " << aWavYR << std::endl;
00381 #endif
00382 
00383         // store descriptor
00384         ioKeyPoint._vec[j++] = aWavXR;
00385         ioKeyPoint._vec[j++] = aWavYR;
00386         /*
00387         if (aWavXR > 0) {
00388                 ioKeyPoint._vec[j++] = aWavXR;
00389                 ioKeyPoint._vec[j++] = 0;
00390         } else {
00391                 ioKeyPoint._vec[j++] = 0;
00392                 ioKeyPoint._vec[j++] = -aWavXR;
00393         }
00394         if (aWavYR > 0) {
00395                 ioKeyPoint._vec[j++] = aWavYR;
00396                 ioKeyPoint._vec[j++] = 0;
00397         } else {
00398                 ioKeyPoint._vec[j++] = 0;
00399                 ioKeyPoint._vec[j++] = -aWavYR;
00400         }
00401         */
00402         if (i != 0)
00403         {
00404             ioKeyPoint._vec[j++] = meanGray - middleMean;
00405         }
00406     }
00407 #ifdef DEBUG_DESC
00408     dlog.close();
00409 #endif
00410 
00411 }

Generated on 29 Aug 2015 for Hugintrunk by  doxygen 1.4.7