Descriptor.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (C) 2007 by Zoran Mesec   *
00003  *   zoran.mesec@gmail.com   *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  ***************************************************************************/
00020 
00021 #include <stdio.h>
00022 #include <iostream>
00023 #include <fstream>
00024 #include <stdlib.h>
00025 #include "Descriptor.h"
00026 
00027 using namespace std;
00028 
00029 Descriptor::Descriptor(APImage* i,HessianDetector* hessianDetector) {
00030     this->image=i;
00031     this->hd=hessianDetector;
00032 }
00033 
00034 void Descriptor::setPoints(vector<vector<int> >* pts) {
00035     this->interestPoints=pts;
00036 }
00037 
00038 // dangelo: this does not compile with MSVC and seems to be unused anyway...
00039 #if 0
00040 void Descriptor::orientate() {
00041     vector<vector<int> >::iterator iter1 = this->interestPoints->begin();
00042     double regionSize;
00043     int kernelSize;
00044 
00045     //int angle;
00046     int max;
00047     int orientation;
00048     double step;
00049 
00050     //coordinates of the sample
00051     int pointX;
00052     int pointY;
00053 
00054 /*
00055 S0: if dx>0 and dy>0 and |dy|>|dx|
00056 S1: if dx>0 and dy>0 and |dy|<=|dx|
00057 S2: if dx>0 and dy<=0 and |dy|<=|dx|
00058 S3: if dx>0 and dy<=0 and |dy|>|dx|
00059 S4: if dx<=0 and dy<=0 and |dy|>|dx|
00060 S5: if dx<=0 and dy<=0 and |dy|<=|dx|
00061 S6: if dx<=0 and dy>0 and |dy|<=|dx|
00062 S7: if dx<=0 and dy>0 and |dy|>|dx|
00063  */
00064     double Sx[NR_ANGLE_BINS];
00065     double Sy[NR_ANGLE_BINS];
00066 
00067 
00068 
00069     double dx;
00070     double dy;
00071 
00072     int pointCount=0;
00073 
00074      while( iter1 != this->interestPoints->end()) { //loop over every interest point
00075          vector<int > interestPoint=*iter1;
00076         //cout << interestPoint[0]<<","<<interestPoint[1]<<":";
00077         regionSize=_getMaxima(interestPoint[0], interestPoint[1])*3;
00078         kernelSize=vigra::round(_getMaxima(interestPoint[0], interestPoint[1])*2);
00079         step=_getMaxima(interestPoint[0], interestPoint[1]);
00080 
00081         for(int i=0; i<NR_ANGLE_BINS;i++) {
00082             Sx[i]=0;
00083             Sy[i]=0;
00084         }
00085 
00086         /*for(int i=interestPoint[0]-regionSize;i<=interestPoint[0]+regionSize;i++) {
00087             for(int j=interestPoint[1]-regionSize; j<=interestPoint[1]+regionSize;j++) {
00088                 cout << "(" <<this->image->getPixel(i,j)<< "," << this->image->getIntegralPixel(i,j) <<") ";
00089             }
00090             cout << "\n";
00091         }*/
00092 
00093         for(double i=(interestPoint[0]-regionSize);i<=(interestPoint[0]+regionSize);i+=step) {
00094             for(double j=(interestPoint[1]-regionSize); j<=(interestPoint[1]+regionSize);j+=step) {
00095 
00096                 //pixels have integer values
00097                 pointX=vigra::round(i);
00098                 pointY=vigra::round(j);
00099                 //pointX and pointY are the coordinates of the sample point
00100 
00101                 double dist=_euclidianDistance(interestPoint[0],interestPoint[1], pointX,pointY);
00102                 if(dist<=regionSize) {   //if point is in the circle
00103 
00104                     double weight=_gaussWeighting(interestPoint[0]-pointX,interestPoint[1]-pointY,2.5*_getMaxima(pointX,pointY));
00105 
00106                     //x directon
00107                     dx=this->image->getRegionSum(pointX-kernelSize, pointY,pointX+kernelSize,pointY+kernelSize);
00108                     dx-=this->image->getRegionSum(pointX-kernelSize, pointY-kernelSize,pointX+kernelSize,pointY);
00109                     dx*=weight;
00110 
00111                     //y direction
00112                     dy=this->image->getRegionSum(pointX,pointY-kernelSize,pointX+kernelSize, pointY+kernelSize);
00113                     dy-=this->image->getRegionSum(pointX-kernelSize,pointY-kernelSize,pointX, pointY+kernelSize);
00114                     dy*=weight;
00115 
00116                     if(dx>0) {
00117                         if(dy>0) {
00118                             if(fabs(dy)>fabs(dx)) {
00119                                 Sx[0]+=dx;
00120                                 Sy[0]+=dy;
00121                             } else {
00122                                 Sx[1]+=dx;
00123                                 Sy[1]+=dy;
00124                             }
00125                         } else {
00126                             if(fabs(dy)<=fabs(dx)) {
00127                                 Sx[2]+=dx;
00128                                 Sy[2]+=dy;
00129                             } else {
00130                                 Sx[3]+=dx;
00131                                 Sy[3]+=dy;
00132                             }
00133                         }
00134                     } else {
00135                         if(dy<=0) {
00136                             if(fabs(dy)>fabs(dx)) {
00137                                 Sx[4]+=dx;
00138                                 Sy[4]+=dy;
00139                             } else {
00140                                 Sx[5]+=dx;
00141                                 Sy[5]+=dy;
00142                             }
00143                         } else {
00144                             if(fabs(dy)<=fabs(dx)) {
00145                                 Sx[6]+=dx;
00146                                 Sy[6]+=dy;
00147                             } else {
00148                                 Sx[7]+=dx;
00149                                 Sy[7]+=dy;
00150                             }
00151                         }
00152                     }
00153 
00154                 } //if point in circle
00155             } //for j
00156         } //for i
00157         max=0;
00158         orientation=0;
00159 
00160         for(int i=0; i<NR_ANGLE_BINS;i++) {
00161             double tmp = _euclidianDistance(0,0,Sx[i],Sy[i]);
00162             if(tmp>max) {
00163                     max = tmp;
00164                     orientation=i;
00165             }
00166         }
00167 
00168         orientation=vigra::round(atan2(Sy[orientation],Sx[orientation]) * 180 / PI);
00169         //if(orientation<0) orientation=360+orientation;
00170 
00171         this->orientations[pointCount]=orientation;
00172 
00173         this->image->drawCircle(interestPoint[1],interestPoint[0],1);
00174         this->image->drawLine(interestPoint[1],interestPoint[0],interestPoint[1]+vigra::round(sin((double)orientation)*10),interestPoint[0]+vigra::round(cos((double)orientation)*10));
00175         //cout << orientation << "\n";
00176        iter1++;
00177        //pointCount++;
00178      }
00179     this->image->show();
00180 }
00181 
00182 #endif
00183 
00184 double Descriptor::_gaussWeighting(int x, int y, double stdev) {
00185 return (1/(2*PI*pow(stdev,2.0)))*exp(-(pow((double)x,2.0)+pow((double)y,2.0))/(2*stdev));
00186 }
00187 double Descriptor::_getMaxima(int x,int y) {
00188     return this->hd->getMaxima(x,y);
00189     //return 1.6;
00190 }
00191 double Descriptor::_euclidianDistance(int x1, int y1, int x2, int y2) {
00192  return sqrt(pow((double)(x1-x2),2.0)+pow((double)(y1-y2),2.0));
00193 }
00194 void Descriptor::_GaborResponse(int x,int y, double maxima, double* descriptor) {
00195     //
00196     double factor = maxima/HD_INITIAL_SCALE;
00197    /*vector<int> d=*descriptor;
00198     d[0]++;*/
00199 
00207     //relative coordinates(from the center of the kernel) and weights
00208     //for points of local maxima for first gabor filter(x,y,weight)
00209     double gabor11[6][3]={
00210     {-1,0,0.839304},
00211     {1,0,-0.839304},
00212     {-3,0,-0.190826},
00213     {3,0,0.190826},
00214     {-6,0,0.0105653},
00215     {6,0,-0.0105653}
00216     };
00217     //second...orientation=0, wavelength=10
00218     double gabor12[6][3]={
00219     {-2,0,0.839304},
00220     {2,0,-0.839304},
00221     {-7,0,-0.20568},
00222     {7,0,0.20568},
00223     {-11,0,0.0133981},
00224     {11,0,-0.0133981}
00225     };
00226 
00227     double gabor13[10][3]={
00228     {-2,0,0.945959},
00229     {2,0,-0.945959},
00230     {-6,0,-0.606531},
00231     {6,0,0.606531},
00232     {-10,0,0.249352},
00233     {10,0,-0.249352},
00234     {-13,0,0.0676238},
00235     {13,0,-0.0676238},
00236     {-17,0,0.0127725},
00237     {17,0,-0.0127725}
00238     };
00239 
00240     //orientation=PI/4, wavelength=5
00241     double gabor21[6][3]={
00242     {-1,-1,0.762278},
00243     {-3,-2,-0.201919},
00244     {2,3,0.201919},
00245     {1,1,0.762278},
00246     {3,5,0.0134254},
00247     {-5,-3,-0.0134254}
00248     };
00249 
00250     //orientation=PI/4, wavelength=10
00251     double gabor22[6][3]={
00252     {-2,-1,0.844207},
00253     {2,1,-0.844207},
00254     {-5,-4,-0.213173},
00255     {5,4,0.213173},
00256     {-9,-7,0.013459},
00257     {9,7,-0.013459}
00258     };
00259 
00260     double gabor23[10][3]={
00261     {-2,-1,0.935087},
00262     {2,1,-0.935087},
00263     {-4,-4,-0.618035},
00264     {4,4,0.618035},
00265     {-7,-7,0.255577},
00266     {7,7,-0.255577},
00267     {-9,-10,-0.0736175},
00268     {10,9,0.0736175},
00269     {-12,-12,0.0126482},
00270     {12,12,-0.0126482}
00271     };
00272 
00273     double gabor31[6][3]={
00274     {0,-1,0.839304},
00275     {0,1,-0.839304},
00276     {0,-3,-0.190826},
00277     {0,3,0.190826},
00278     {0,6,-0.0105653},
00279     {0,-6,0.0105653}
00280     };
00281     double gabor32[6][3]={
00282     {0,2,-0.839304},
00283     {0,7,0.20568},
00284     {0,11,-0.0133981},
00285     {0,-2,0.839304},
00286     {0,-7,-0.20568},
00287     {0,-11,0.0133981}
00288     };
00289     double gabor33[10][3]={
00290     {0,-2,0.945959},
00291     {0,2,-0.945959},
00292     {0,-6,-0.606531},
00293     {0,6,0.606531},
00294     {0,-10,0.249352},
00295     {0,10,-0.249352},
00296     {0,-13,-0.0676238},
00297     {0,13,0.0676238},
00298     {0,-17,-0.0127725},
00299     {0,17,0.0127725}
00300     };
00301 
00302     double gabor41[6][3]={
00303     {1,-1,0.762278},
00304     {-2,+3,0.201919},
00305     {2,-3,-0.201919},
00306     {-1,1,-0.762278},
00307     {3,-5,0.0134254},
00308     {-3,5,-0.0134254}
00309     };
00310     double gabor42[6][3]={
00311     {-1,2,-0.844207},
00312     {2,-1,0.844207},
00313     {4,-5,-0.213173},
00314     {-5,4,0.213173},
00315     {8,-8,0.013459},
00316     {-8,8,-0.013459}
00317     };
00318     double gabor43[10][3]={
00319     {2,-1,0.935087},
00320     {-1,2,-0.935087},
00321     {4,-4,-0.618035},
00322     {-4,4,0.618035},
00323     {7,-7,0.255577},
00324     {-7,7,-0.255577},
00325     {9,-10,-0.0736175},
00326     {-10,9,0.0736175},
00327     {12,-12,0.0126482},
00328     {-12,12,-0.0126482}
00329     };
00330 
00331     //orientation PI/8
00332     double gabor53[10][3]={
00333     {-2,0,0.946801},
00334     {2,0,-0.946801},
00335     {-5,-3,-0.619484},
00336     {3,5,0.619484},
00337     {-7,-8,0.263348},
00338     {8,7,-0.263348},
00339     {-12,-6,0.0735327},
00340     {6,12,-0.0735327},
00341     {-15,-9,0.0133379},
00342     {9,15,-0.0133379}
00343     };
00344 
00345     //orientation 3*(PI/8)
00346     double gabor63[10][3]={
00347     {0,-2,0.946801},
00348     {0,2,-0.946801},
00349     {-3,-5,-0.619484},
00350     {5,3,0.619484},
00351     {-8,-7,0.263348},
00352     {7,8,-0.263348},
00353     {-6,-12,0.0735327},
00354     {12,6,-0.0735327},
00355     {-9,-15,0.0133379},
00356     {15,9,-0.0133379}
00357     };
00358 
00359     //orientation 5*PI/8
00360     double gabor73[10][3]={
00361     {0,-2,0.946801},
00362     {0,2,-0.946801},
00363     {3,-5,-0.619484},
00364     {-3,5,0.619484},
00365     {8,-7,0.263348},
00366     {-8,8,-0.263348},
00367     {6,-12,0.0735327},
00368     {-6,12,-0.0735327},
00369     {9,-15,0.0133379},
00370     {-9,15,-0.0133379}
00371     };
00372 
00373     //orientation 7*PI/8
00374     double gabor83[10][3]={
00375     {2,0,0.946801},
00376     {-2,0,-0.946801},
00377     {5,-3,-0.619484},
00378     {-5,3,0.619484},
00379     {7,-8,0.263348},
00380     {-7,8,-0.263348},
00381     {12,-6,0.0735327},
00382     {-12,6,-0.0735327},
00383     {15,-9,0.0133379},
00384     {-15,9,-0.0133379}
00385     };
00386 
00387     //stdev=6, wavelength=4, phase=0, orientation=0
00388     double gabor91[17][3]={
00389     {-16,0, 0.0285655},
00390     {-14,0, -0.0657285},
00391     {-12,0, 0.135335},
00392     {-10,0, -0.249352},
00393     {-8,0, 0.411112},
00394     {-6,0, -0.606531},
00395     {-4,0, 0.800737},
00396     {-2,0, -0.945959},
00397     {0,0, 1},
00398     {16,0, 0.0285655},
00399     {14,0, -0.0657285},
00400     {12,0, 0.135335},
00401     {10,0, -0.249352},
00402     {8,0, 0.411112},
00403     {6,0, -0.606531},
00404     {4,0, 0.800737},
00405     {2,0, -0.945959}
00406     };
00407 
00408     //stdev=6, wavelength=4, phase=0, orientation=PI/8
00409     double gabor101[17][3]={
00410     {-16,-3,0.0292416},
00411     {-14,-2, -0.0656082},
00412     {-12,-2, 0.138167},
00413     {-10,-2, -0.248923},
00414     {-8,-1, 0.404746},
00415     {-6,-1, -0.609707},
00416     {-4,-1, 0.787721},
00417     {-2,0, -0.926472},
00418     {0,0, 1},
00419     {16,3,0.0292416},
00420     {14,2, -0.0656082},
00421     {12,2, 0.138167},
00422     {10,2, -0.248923},
00423     {8,1, 0.404746},
00424     {6,1, -0.609707},
00425     {4,1, 0.787721},
00426     {2,0, -0.926472}
00427     };
00428     //stdev=6, wavelength=4, phase=0, orientation=PI/4
00429     double gabor111[17][3]={
00430     {-11,-11, 0.026607},
00431     {-10,-10, -0.0606333},
00432     {-9,-8, 0.134318},
00433     {-7,-7, -0.253187},
00434     {-6,-5, 0.405626},
00435     {-4,-4, -0.550271},
00436     {-3,-3, 0.722915},
00437     {-2,-1, -0.922342},
00438     {0,0, 1},
00439     {11,11, 0.026607},
00440     {10,10, -0.0606333},
00441     {8,9, 0.134318},
00442     {7,7, -0.253187},
00443     {5,6, 0.405626},
00444     {4,4, -0.550271},
00445     {3,3, 0.722915},
00446     {1,2, -0.922342}
00447     };
00448 
00449     //3PI/8
00450      double gabor121[17][3]={
00451     {-5,-15, 0.0295778},
00452     {-5,-13, -0.0672137},
00453     {-2,-12, 0.138167},
00454     {-4,-9, -0.252577},
00455     {-1,-8, 0.404746},
00456     {-3,-5, -0.588397},
00457     {-1,-4, 0.787721},
00458     {-0,-2, -0.926472},
00459     {0,0, 1},
00460     {5,15, 0.0295778},
00461     {5,13, -0.0672137},
00462     {2,12, 0.138167},
00463     {4,9, -0.252577},
00464     {1,8, 0.404746},
00465     {3,5, -0.588397},
00466     {1,4, 0.787721},
00467     {0,2, -0.926472}
00468     };
00469 
00470     double gabor131[17][3]={
00471     {0,-16, 0.0285655},
00472     {0,-14, -0.0657285},
00473     {0,-12, 0.135335},
00474     {0,-10, -0.249352},
00475     {0,-8, 0.411112},
00476     {0,-6, -0.606531},
00477     {0,-4, 0.800737},
00478     {0,-2, -0.945959},
00479     {0,0, 1},
00480     {0,16, 0.0285655},
00481     {0,14, -0.0657285},
00482     {0,12, 0.135335},
00483     {0,10, -0.249352},
00484     {0,8, 0.411112},
00485     {0,6, -0.606531},
00486     {0,4, 0.800737},
00487     {0,2, -0.945959}
00488     };
00489 
00490     //5PI/8
00491     double gabor141[17][3]={
00492     {10,-13, 0.0296244},
00493     {5,-14, -0.0672137},
00494     {7,-10, 0.137794},
00495     {4,-9, -0.252577},
00496     {6,-6, 0.411811},
00497     {1,-6, -0.606531},
00498     {3,-3, 0.801129},
00499     {3,-1, -0.93537},
00500     {0,0, 1},
00501     {-10,13, 0.0296244},
00502     {-5,14, -0.0672137},
00503     {-7,10, 0.137794},
00504     {-4,9, -0.252577},
00505     {-6,6, 0.411811},
00506     {-1,6, -0.606531},
00507     {-3,-3, 0.801129},
00508     {-3,1, -0.93537}
00509     };
00510 
00511     //6PI/8
00512     double gabor151[17][3]={
00513     {11,-11, 0.026607},
00514     {10,-10, -0.0606333},
00515     {9,-8, 0.134318},
00516     {7,-7, -0.253187},
00517     {6,-5, 0.405626},
00518     {4,-4, -0.550271},
00519     {3,-3, 0.722915},
00520     {2,-1, -0.922342},
00521     {0,0, 1},
00522     {-11,11, 0.026607},
00523     {-10,10, -0.0606333},
00524     {-8,9, 0.134318},
00525     {-7,7, -0.253187},
00526     {-5,6, 0.405626},
00527     {-4,4, -0.550271},
00528     {-3,3, 0.722915},
00529     {-1,2, -0.922342}
00530     };
00531 
00532     double gabor161[17][3] = {
00533     {15,-5, 0.0295778},
00534     {13,-5, -0.0672137},
00535     {12,-2, 0.138167},
00536     {9,-4, -0.252577},
00537     {8,-1, 0.404746},
00538     {5,-3, -0.588397},
00539     {4,-1, 0.787721},
00540     {0,-2, -0.926472},
00541     {0,0, 1},
00542     {-15,5, 0.0295778},
00543     {-13,5, -0.0672137},
00544     {-12,2, 0.138167},
00545     {-9,4, -0.252577},
00546     {-8,1, 0.404746},
00547     {-5,3, -0.588397},
00548     {-4,1, 0.787721},
00549     {-2,0, -0.926472}
00550     };
00551 
00552 
00553     for(int i=0; i<6; i++) {
00554 
00555         descriptor[0]+=this->image->getPixel(vigra::round(x+factor*gabor11[i][0]),vigra::round(y+factor*gabor11[i][1]))*gabor11[i][2];
00556         descriptor[2]+=this->image->getPixel(vigra::round(x+factor*gabor21[i][0]),vigra::round(y+factor*gabor21[i][1]))*gabor21[i][2];
00557         descriptor[4]+=this->image->getPixel(vigra::round(x+factor*gabor31[i][0]),vigra::round(y+factor*gabor31[i][1]))*gabor31[i][2];
00558         descriptor[6]+=this->image->getPixel(vigra::round(x+factor*gabor41[i][0]),vigra::round(y+factor*gabor41[i][1]))*gabor41[i][2];
00559 
00560         descriptor[1]+=this->image->getPixel(vigra::round(x+factor*gabor12[i][0]),vigra::round(x+factor*gabor12[i][1]))*gabor12[i][2];
00561         descriptor[3]+=this->image->getPixel(vigra::round(x+factor*gabor22[i][0]),vigra::round(x+factor*gabor22[i][1]))*gabor22[i][2];
00562         descriptor[5]+=this->image->getPixel(vigra::round(x+factor*gabor32[i][0]),vigra::round(x+factor*gabor32[i][1]))*gabor32[i][2];
00563         descriptor[7]+=this->image->getPixel(vigra::round(x+factor*gabor42[i][0]),vigra::round(x+factor*gabor42[i][1]))*gabor42[i][2];
00564 
00565     }
00566     for(int i=0; i<10;i++) {
00567 
00568         descriptor[8]+=this->image->getPixel(vigra::round(x+factor*gabor13[i][0]),vigra::round(y+factor*gabor13[i][1]))*gabor13[i][2];
00569         descriptor[9]+=this->image->getPixel(vigra::round(x+factor*gabor23[i][0]),vigra::round(y+factor*gabor23[i][1]))*gabor23[i][2];
00570         descriptor[10]+=this->image->getPixel(vigra::round(x+factor*gabor33[i][0]),vigra::round(y+factor*gabor33[i][1]))*gabor33[i][2];
00571         descriptor[11]+=this->image->getPixel(vigra::round(x+factor*gabor43[i][0]),vigra::round(y+factor*gabor43[i][1]))*gabor43[i][2];
00572 
00573         descriptor[12]+=this->image->getPixel(vigra::round(x+factor*gabor53[i][0]),vigra::round(x+factor*gabor53[i][1]))*gabor53[i][2];
00574         descriptor[13]+=this->image->getPixel(vigra::round(x+factor*gabor63[i][0]),vigra::round(x+factor*gabor63[i][1]))*gabor63[i][2];
00575         descriptor[14]+=this->image->getPixel(vigra::round(x+factor*gabor73[i][0]),vigra::round(x+factor*gabor73[i][1]))*gabor73[i][2];
00576         descriptor[15]+=this->image->getPixel(vigra::round(x+factor*gabor83[i][0]),vigra::round(x+factor*gabor83[i][1]))*gabor83[i][2];
00577     }
00578     for(int i=0; i<17;i++) {
00579 
00580         descriptor[16]+=this->image->getPixel(vigra::round(x+factor*gabor91[i][0]),vigra::round(y+factor*gabor91[i][1]))*gabor91[i][2];
00581         descriptor[17]+=this->image->getPixel(vigra::round(x+factor*gabor101[i][0]),vigra::round(y+factor*gabor101[i][1]))*gabor101[i][2];
00582         descriptor[18]+=this->image->getPixel(vigra::round(x+factor*gabor111[i][0]),vigra::round(y+factor*gabor111[i][1]))*gabor111[i][2];
00583         descriptor[19]+=this->image->getPixel(vigra::round(x+factor*gabor121[i][0]),vigra::round(y+factor*gabor121[i][1]))*gabor121[i][2];
00584 
00585         descriptor[20]+=this->image->getPixel(vigra::round(x+factor*gabor13[i][0]),vigra::round(y+factor*gabor131[i][1]))*gabor131[i][2];
00586         descriptor[21]+=this->image->getPixel(vigra::round(x+factor*gabor141[i][0]),vigra::round(y+factor*gabor141[i][1]))*gabor141[i][2];
00587         descriptor[22]+=this->image->getPixel(vigra::round(x+factor*gabor151[i][0]),vigra::round(y+factor*gabor151[i][1]))*gabor151[i][2];
00588         descriptor[23]+=this->image->getPixel(vigra::round(x+factor*gabor161[i][0]),vigra::round(y+factor*gabor161[i][1]))*gabor161[i][2];
00589 
00590     }
00591 }
00592 
00593 /*void Descriptor::_RegionResponse(int x,int y, double maxima, double* descriptor) {
00594     int add;
00595     int sum1;
00596     int sum2;
00597     int sum3;
00598     int sum4;
00599     int cnt=0;
00600 
00601     for(int i=2; i<=32;i*=2) {
00602         add=vigra::round(maxima*i);
00603         sum1=this->image->getRegionSum(x-add,y,x,y+add);
00604         sum2=-1*this->image->getRegionSum(x,y,x+add,y+add);
00605         sum3=this->image->getRegionSum(x,y-add,x+add,y);
00606         sum4=-1*this->image->getRegionSum(x-add,y-add,x,y);
00607         descriptor[cnt+0]+=(sum1+sum3)/pow((double)i,4.0);
00608         descriptor[cnt+1]+=(sum2+sum4)/pow(i,4);
00609         descriptor[cnt+2]+=(sum1+sum2+sum3+sum4)/pow(i,4);
00610         descriptor[cnt+3]+=abs(sum1+sum2+sum3+sum4)/pow(i,4);
00611         cnt+=5;
00612     }
00613 }*/
00614 
00615 void Descriptor::_ShapeContext(int x,int y, double maxima, double* descriptor) {
00616 
00617 
00618     /*for(int i=2; i<=32;i*=2) {
00619         add=vigra::round(maxima*i);
00620         sum1=this->image->getRegionSum(x-add,y,x,y+add);
00621         sum2=-1*this->image->getRegionSum(x,y,x+add,y+add);
00622         sum3=this->image->getRegionSum(x,y-add,x+add,y);
00623         sum4=-1*this->image->getRegionSum(x-add,y-add,x,y);
00624         descriptor[cnt+0]+=(sum1+sum3)/pow((double)i,4.0);
00625         descriptor[cnt+1]+=(sum2+sum4)/pow(i,4);
00626         descriptor[cnt+2]+=(sum1+sum2+sum3+sum4)/pow(i,4);
00627         descriptor[cnt+3]+=abs(sum1+sum2+sum3+sum4)/pow(i,4);
00628         cnt+=5;
00629     }*/
00630 }
00631 
00632 void Descriptor::createDescriptors() {
00633 
00634 // TODO (zoran#8#): Current version detects edges using canny algorithm in the region around the interest point -  here is a lot of redundancy. \
00635 // Need to change: using vigras cannyEdgeImage first run edge detection on the entire picture and then create desciptors.
00636     cout << "Creating descriptors..."<< endl;
00637     int regionSize;
00638     int i=0;
00639     double orientation;
00640     this->interestPoints->resize(this->interestPoints->size());
00641     cout << this->interestPoints->size()<<endl;
00642     vector<vector<int> >::iterator iter1 = this->interestPoints->begin();
00643          while( iter1 != this->interestPoints->end() ) { //loop over every interest point
00644             vector<int > interestPoint=*iter1;
00645             i++;
00646             vector<double> descriptor;
00647             for(int i=0;i<DESCRIPTOR_SIZE;i++) { //initialize array
00648                 descriptor.push_back(0);
00649             }
00650 
00651             regionSize=vigra::round(15*this->_getMaxima(interestPoint[0], interestPoint[1]));
00652             double b0,b1;
00653             b0=6*this->_getMaxima(interestPoint[0], interestPoint[1]);
00654             b1=11*this->_getMaxima(interestPoint[0], interestPoint[1]);
00655 
00656             int xStart=(interestPoint[0]-regionSize);
00657             int xEnd=(interestPoint[0] +regionSize);
00658             int yStart=(interestPoint[1]-regionSize);
00659             int yEnd=(interestPoint[1]+regionSize);
00660 
00661             // only analyse pixels in image
00662             int xStartImg = xStart < 0 ? 0 : xStart;
00663             int xEndImg = xEnd >= this->image->getWidthBW() ? this->image->getWidthBW() -1  : xEnd;
00664             int yStartImg = yStart < 0 ? 0 : yStart;
00665             int yEndImg = yEnd >= this->image->getHeightBW() ? this->image->getHeightBW() -1 : yEnd;
00666 
00667 /*
00668                             double result;
00669                             result = atan (gy/gx) * 180 / PI;
00670                             result+=90;
00671 
00672                             if(result<22.5 && result>=157.5) cp->direction=0;
00673                             else if(result>=22.5 && result<67.5) cp->direction=1;
00674                             else if(result>=67.5 && result<112.5) cp->direction=2;
00675                             else if(result>=112.5 && result<157.5) cp->direction=3;
00676 
00677 */
00678             // empty edgel list
00679             std::vector<vigra::Edgel> edgels;
00680 
00681             //cout << "Maxima:" << this->_getMaxima(interestPoint[0], interestPoint[1]) << endl;
00682             //cout << "Point:" << interestPoint[0] << "," << interestPoint[1]<< endl;
00683 
00684             //
00685             vector<int > iPointRelative = interestPoint;
00686             iPointRelative[0]=regionSize;
00687             iPointRelative[1]=regionSize;
00688 
00689             // find edgels at scale of the interest point
00690             vigra::cannyEdgelList(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStartImg, yStartImg),
00691                                                 this->image->imgBW->upperLeft() + vigra::Diff2D(xEndImg, yEndImg)),
00692                                   edgels, this->_getMaxima(interestPoint[0], interestPoint[1]),iPointRelative);
00693             /*this->image->_cannyEdgelList1(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStart, yStart),
00694                                                 this->image->imgBW->upperLeft() + vigra::Diff2D(xEnd, yEnd)),
00695                                   edgels, this->_getMaxima(interestPoint[0], interestPoint[1]), &interestPoint);*/
00696             //cout << "Size:" << edgels.size() << endl;
00697 
00698             //first edgel element holds orientation assignment for interest point
00699             vector<vigra::Edgel>::iterator iter2 = edgels.begin();
00700             vigra::Edgel edgePoint=*iter2;
00701             iter2++;
00702             //cout << "Orientation:" << edgePoint.orientation << ", strenght:" << edgePoint.strength << endl;
00703             orientation = edgePoint.orientation;
00704             if(orientation<0) orientation+=2*PI;
00705 
00706             double a0,a1,a2,a3;
00707             a0=fmod(0+orientation,2*PI);
00708             a1=fmod(PI/2+orientation,2*PI);
00709             a2=fmod(PI+orientation,2*PI);
00710             a3=fmod((3*PI)/2+orientation,2*PI);
00711 
00712             //int count =0;
00713              while( iter2 != edgels.end()) { //loop over every canny pixel
00714                 edgePoint=*iter2;
00715 
00716                 //TODO discard edgels that are not in circle
00717                 double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x+xStartImg, edgePoint.y+yStartImg);
00718 
00719                 if(distance>regionSize) {
00720                     iter2++; continue;
00721                 }
00722                 //count++;
00723                 //cout << edgePoint.x << ","<<edgePoint.y << "; orientation:"<<edgePoint.orientation<<",strength:"<<edgePoint.strength<<endl;
00724                 int sector;
00725                 if(edgePoint.orientation<0) edgePoint.orientation+=2*PI;
00726 
00727                 int orientSector;
00728 
00729                 //orientation is quantized into 4 bins
00730                 if(edgePoint.orientation>=0 && edgePoint.orientation<(PI/2)) orientSector=0;
00731                 else if(edgePoint.orientation>=(PI/2) && edgePoint.orientation<PI) orientSector=1;
00732                 else if(edgePoint.orientation>=PI && edgePoint.orientation<(3*PI/2)) orientSector=2;
00733                 else if(edgePoint.orientation>=(3*PI/2)) orientSector=3;
00734                 else cerr<< "Orientation1:"<<edgePoint.orientation<<endl;
00735 
00736                 double result = atan2(regionSize-edgePoint.y,edgePoint.x-regionSize);
00737                 if(result<0.0) result+=2*PI;
00738 
00739                 //location is quantized into 9 bins
00740                 if(result>=a0 && result<a1) sector=0;
00741                 else if(result>=a1 && result<a2) sector=1;
00742                 else if(result>=a2 && result<a3) sector=2;
00743                 else if(result>=a3 || result < a0) sector=3;
00744                 else cerr<< "Orientation2:"<<edgePoint.orientation<<endl;
00745 
00746                 if(distance<=b0) descriptor[orientSector*9 + 8]+=edgePoint.strength;
00747                 else if(distance>b0 && distance<=b1) descriptor[orientSector*9 + sector*2+0]+=edgePoint.strength;
00748                 else if(distance>b1) descriptor[orientSector*9 + sector*2+1]+=edgePoint.strength;
00749 
00750                 /*double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x,edgePoint.y);
00751                 if(distance<=27) descriptor[orientSector*9 + 8]+=edgePoint.strength;
00752                 else if(distance>27 && distance<=50) descriptor[orientSector*9 + sector*2+0]+=edgePoint.strength;
00753                 else if(distance>50) descriptor[orientSector*9 + sector*2+1]+=edgePoint.strength;
00754 */
00755 
00756                 iter2++;
00757              }
00758                 //cout << "Size: "<<edgels.size()<< "; count: "<<count<<endl;
00759             this->descriptors.push_back(descriptor);
00760             iter1++;
00761          }
00762 }
00763 vector<vector<double> >* Descriptor::getDescriptors() {
00764     return &this->descriptors;
00765 }
00766 bool Descriptor::printDescriptors(string name)
00767 {
00768     ofstream o (name.c_str());
00769         if (o.is_open()) {
00770         o << DESCRIPTOR_SIZE << endl;
00771         o << interestPoints->size() << endl;
00772         vector<vector<int> >::iterator iter1 = interestPoints->begin();
00773         vector<vector<double> >::iterator iterDesc = descriptors.begin();
00774 
00775         //int c=0;
00776         while( iter1 != interestPoints->end()) {
00777              vector<int > tmp2=*iter1;
00778              double r = _getMaxima(tmp2[0], tmp2[1])*HD_INIT_KERNEL_SIZE;
00779              o <<tmp2[1]<<" "<<tmp2[0]<<" " << 1/(r*r) << " " << 0 << " " <<  1/(r*r) << " ";
00780              for (vector<double>::iterator it = iterDesc->begin(); it != iterDesc->end(); ++it) {
00781                 o << *it << " ";
00782              }
00783              o << endl;
00784              iter1++;
00785              iterDesc++;
00786         }
00787         o.close();
00788         return true;
00789     } else return false;
00790 }
00791 
00792 bool Descriptor::generateAutopanoXML(string name) {
00793     cout << "Generating XML file: "<<name<<endl;
00794         ofstream xmlFile (name.c_str());
00795         if (xmlFile.is_open())
00796         {
00797     xmlFile << "<?xml version=\"1.0\" encoding=\"utf-8\"?>"<< endl;
00798     xmlFile << "<KeypointXMLList xmlns:xsd=\"http://www.w3.org/2001/XMLSchema\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\">\n";
00799     xmlFile << "  <XDim>"<< this->image->getWidth() <<"</XDim>"<< endl;
00800     xmlFile << "  <YDim>"<< this->image->getHeight() <<"</YDim>"<< endl;
00801     xmlFile << "  <ImageFile>"<< this->image->getPath() <<"</ImageFile>"<< endl;
00802     xmlFile << "  <Arr>"<< endl;
00803 
00804     vector<vector<int> >::iterator iter1 = interestPoints->begin();
00805     vector<vector<double> >::iterator iterDesc = descriptors.begin();
00806 
00807     int c=0;
00808     while( iter1 != interestPoints->end()) {
00809          vector<int > tmp2=*iter1;
00810 
00811          xmlFile << "    <KeypointN>"<< endl;
00812          xmlFile << "      <X>"<< tmp2[0] <<"</X>"<<endl;
00813          xmlFile << "      <Y>"<< tmp2[1] <<"</Y>"<<endl;
00814          xmlFile << "      <Scale>"<< this->_getMaxima(tmp2[0],tmp2[1]) <<"</Scale>"<<endl;
00815          xmlFile << "      <Orientation>"<< 0 <<"</Orientation>"<<endl;
00816          xmlFile << "      <Dim>"<< DESCRIPTOR_SIZE <<"</Dim>"<<endl;
00817          xmlFile << "      <Descriptor>"<<endl;
00818 
00819          for (vector<double>::iterator it = iterDesc->begin(); it != iterDesc->end(); ++it) {
00820                  xmlFile << "        <int>"<< (int)*it << "</int>" << endl;
00821          }
00822          xmlFile << "      </Descriptor>"<<endl;
00823          xmlFile << "    </KeypointN>"<< endl;
00824 
00825          iter1++;
00826          iterDesc++;
00827          c++;
00828     }
00829 
00830     xmlFile << "  </Arr>"<< endl;
00831     xmlFile << "</KeypointXMLList>"<< endl;
00832 
00833     xmlFile.close();
00834     return true;
00835   }
00836   else return false;
00837 }
00838 
00839 

Generated on Mon Sep 9 01:25:42 2013 for Hugintrunk by  doxygen 1.3.9.1