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

Generated on Thu Aug 21 01:25:41 2014 for Hugintrunk by  doxygen 1.3.9.1