HessianDetector.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 #include <stdio.h>
00021 #include <iostream>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include <limits>
00025 #include "HessianDetector.h"
00026 
00027 using namespace std;
00028 
00029 HessianDetector::HessianDetector(APImage* i, int nrPoints, CONVOLUTION_TYPE type, int octaves) {
00030     if(octaves>HD_MAX_OCTAVES) return;  // TODO (zoran#1#): Throw error or log status.
00031 
00032         this->image=i;
00033         this->nrPoints=nrPoints;
00034         this->nrOctaves=octaves;
00035         this->convolutionType = type;
00036 
00037     vector<int> point;
00038     point.resize(2);
00039 
00040     this->orderedList.reserve(sizeof(point)*nrPoints);
00041 
00042     //initialize the vectors
00043         determinants.resize(this->image->getWidthBW());
00044         maximas.resize(this->image->getWidthBW());
00045 
00046     const int iMax=std::numeric_limits<int>::max();
00047         for(int i=0;i<this->image->getWidthBW();i++) {
00048                 maximas[i].resize(this->image->getHeightBW());
00049                 determinants[i].resize(this->image->getHeightBW());
00050                 for(int j=0; j<this->image->getHeightBW();j++) {
00051                         maximas[i][j] =0 ;
00052                         determinants[i][j]=iMax;
00053                 }
00054         }
00055 }
00056 
00057 bool HessianDetector::detect() {
00058     if(this->convolutionType==HD_BOX_FILTERS) {
00059      return this->_boxFilterDetect();
00060     } else if(this->convolutionType==HD_SLIDING_WINDOW) {
00061      return this->_slidingWDetect();
00062     }
00063     return false;
00064 }
00065 
00066 bool HessianDetector::_slidingWDetect() {
00067         int height = this->image->getHeightBW();
00068         int width = this->image->getWidthBW();
00069 
00070         int kernelX=9;
00071         int kernelY=9;
00072         int kernelSize=9;
00073         double scale=1.2;
00074 
00075         int gxx[][9]={
00076         1, 7, 29, 68, 90, 68, 29, 7, 1,   //9x9 gauss kernels
00077     4, 26, 106, 247, 327, 247, 106, 26, 4,
00078     5, 33, 133, 310, 411, 310, 133, 33, 5,
00079     -4, -27, -109, -252, -335, -252, -109, -27, -4,
00080     -11, -81, -329, -765, -1013, -765, -329, -81, -11,
00081     -4, -27, -109, -252, -335, -252, -109, -27, -4,
00082     5, 33, 133, 310, 411, 310, 133, 33, 5,
00083     4, 26, 106, 247, 327, 247, 106, 26, 4,
00084     1, 7, 29, 68, 90, 68, 29, 7, 1
00085         };
00086 
00087         int gyy[][9]={
00088     1, 4, 5, -4, -11, -4, 5, 4, 1,
00089     7, 26, 33, -27, -81, -27, 33, 26, 7,
00090     29, 106, 133, -109, -329, -109, 133, 106, 29,
00091     68, 247, 310, -252, -765, -252, 310, 247, 68,
00092     90, 327, 411, -335, -1013, -335, 411, 327, 90,
00093     68, 247, 310, -252, -765, -252, 310, 247, 68,
00094     29, 106, 133, -109, -329, -109, 133, 106, 29,
00095     7, 26, 33, -27, -81, -27, 33, 26, 7,
00096     1, 4, 5, -4, -11, -4, 5, 4, 1
00097     };
00098 
00099     int gxy[][9]={
00100     1, 5, 15, 17, 0, -17, -15, -5, -1,
00101     5, 29, 78, 91, 0, -91, -78, -29, -5,
00102     15, 78, 214, 248, 0, -248, -214, -78, -15,
00103     17, 91, 248, 289, 0, -289, -248, -91, -17,
00104     0, 0, 0, 0, 0, 0, 0, 0, 0,
00105     -17, -91, -248, -289, 0, 289, 248, 91, 17,
00106     -15, -78, -214, -248, 0, 248, 214, 78, 15,
00107     -5, -29, -78, -91, 0, 91, 78, 29, 5,
00108     -1, -5, -15, -17, 0, 17, 15, 5, 1
00109     };
00110 
00111     int offsetX=floor((double)kernelX/2);
00112     int offsetY=floor((double)kernelY/2);
00113     int xStart=0;
00114     int yStart=0;
00115     int pixelSumYY;
00116     int pixelSumXX;
00117     int pixelSumXY;
00118 
00119     int tmpX=0;
00120     int tmpY=0;
00121     int val;
00122     int mappingX;
00123     int mappingY;
00124 
00125         int det;
00126 
00127 
00128     for(int s=0; s<6;s+=1.0) {   //scale-space loop
00129     //cout << s << "\n";
00130     height=image->getHeightBW();
00131     width=image->getWidthBW();
00132     /*cout << height << "," << width << " \n";
00133     cout << " ";*/
00134     for(int i=0;i<height;i++) {
00135         for(int j=0; j<width;j++) {
00136             pixelSumYY=0;
00137             pixelSumXX=0;
00138             pixelSumXY=0;
00139 
00140             xStart=i-offsetX;
00141             yStart=j-offsetY;
00142 
00143             for(int k=0;k<kernelX;k++) {
00144                 tmpX=xStart+k;
00145                 if(tmpX<0) tmpX=height+tmpX;
00146                 if(tmpX>=height) tmpX=tmpX-height;
00147 
00148                 for(int l=0;l<kernelY;l++) {
00149                     tmpY=yStart+l;
00150 
00151                     if(tmpY<0) tmpY=width+tmpY;
00152                     if(tmpY>=height) tmpY=tmpY-height;
00153 
00154                     val=image->getPixel(tmpX,tmpY);
00155 
00156                     pixelSumYY+=gyy[k][l]*val;
00157                     pixelSumXX+=gxx[k][l]*val;
00158                     pixelSumXY+=gxy[k][l]*val;
00159                 }
00160             }
00161 
00162             //determinant of the Hessian Matrix
00163                         det = pixelSumXX*pixelSumYY-pow((double)pixelSumXY,2.0);
00164 
00165             //map coordinates on the scaled image onto coordinates on the original image
00166             mappingX=vigra::round(i*pow(scale,s));
00167             mappingY=vigra::round(j*pow(scale,s));
00168 
00169                         if(det>determinants[mappingX][mappingY]) {
00170                 //if(det>max) max=det;
00171                 determinants[mappingX][mappingY]=det;
00172                 maximas[mappingX][mappingY]=kernelSize*pow(scale,s);
00173                         }
00174         }
00175        //cout << "\n";
00176     }
00177     image->scale(scale);
00178     }   //scale-space
00179 
00180     int xEnd=0;
00181     int yEnd=0;
00182 
00183     for(int i=0;i<image->getHeight();i++) {
00184         if(i>=(image->getHeight()-1)) {
00185             xStart=i-1;
00186             xEnd=0;
00187         } else if(i==0) {
00188             xStart=image->getHeight()-1;
00189             xEnd=i+1;
00190         } else {
00191             xStart=i-1;
00192             xEnd=i+1;
00193         }
00194         for(int j=0; j<image->getWidth();j++) {
00195             if(j>=(image->getWidth()-1)) {
00196                 yStart=j-1;
00197                 yEnd=0;
00198             } else if(j==0) {
00199                 yStart=image->getWidth()-1;
00200                 yEnd=j+1;
00201             } else {
00202                 yStart=j-1;
00203                 yEnd=j+1;
00204             }
00212             if(determinants[i][j]>determinants[xStart][yStart] &&
00213             determinants[i][j]>determinants[xStart][j] &&
00214             determinants[i][j]>determinants[xStart][yEnd] &&
00215             determinants[i][j]>determinants[i][yStart] &&
00216             determinants[i][j]>determinants[i][yEnd] &&
00217             determinants[i][j]>determinants[xEnd][yStart] &&
00218             determinants[i][j]>determinants[xEnd][j] &&
00219             determinants[i][j]>determinants[xEnd][yEnd] ) {
00220                 this->_insertToList(&i,&j);
00221             }
00222 
00223         }
00224         //cout << "\n";
00225     }
00226     return true;
00227 }
00228 
00229 bool HessianDetector::_boxFilterDetect() {
00230 
00231         int height = this->image->getHeightBW();
00232         int width = this->image->getWidthBW();
00233 
00234     cout << "Height:"<<height << ",width:"<<width<<endl;
00235     //cout << "Max vector size:"<< orderedList.max_size() << ", capacity:" << orderedList.capacity() << "\n";
00236 
00237 
00238     const int iMax=std::numeric_limits<int>::max();
00239 
00240     int offsetX=width*0.05;
00241     int offsetY=height*0.05;
00242 
00243 
00244     int xStart=offsetX;
00245     int yStart=offsetY;
00246     int xEnd=width-offsetX;
00247     int yEnd=height-offsetY;
00248 
00249     for(int j=yStart;j<yEnd;j++) {
00250         for(int i=xStart; i<xEnd;i++) {
00251             //_calculateMaxDet(x_coord,y_coord)
00252             _calculateMaxDet(i,j);  //calculate determinant
00253 
00254         }
00255     }
00256 
00257  vector<vector<int> >::iterator it;
00258 
00259   it = orderedList.begin();
00260   double sum=0;
00261   double count=0;
00262     //point (x,y) => (j,i)
00263     for(int j=yStart;j<yEnd;j++) {
00264         for(int i=xStart; i<xEnd;i++) {
00273             if(determinants[i][j]>determinants[i-1][j-1] &&
00274             determinants[i][j]>determinants[i-1][j] &&
00275             determinants[i][j]>determinants[i-1][j+1] &&
00276             determinants[i][j]>determinants[i][j-1] &&
00277             determinants[i][j]>determinants[i][j+1] &&
00278             determinants[i][j]>determinants[i+1][j-1] &&
00279             determinants[i][j]>determinants[i+1][j] &&
00280             determinants[i][j]>determinants[i+1][j+1] &&
00281             determinants[i][j]>300000) {
00282             //this->_insertToList(&i,&j);
00283             //cout << maximas[i][j]<<","<<determinants[i][j]<< ";"<<endl;
00284             /*cout << i<<","<<j<< "\n";*/
00285             vector<int> point;
00286             point.resize(2);
00287 
00288             point[0]=i;
00289             point[1]=j;
00290 
00291             this->orderedList.push_back(point);
00292 
00293             count++;
00294             sum+=determinants[i][j];
00295             //it = orderedList.insert ( it , point);
00296 
00297             }
00298                 //cout << maximas[i][j]<<","<<determinants[i][j]<< " ";
00299         }
00300        //cout << "\n";
00301     }
00302     cout << "Detected points:"<< this->orderedList.size() <<" ";
00303     cout << "Average:"<< sum/count <<endl;
00304     this->_cutPointList(sum/count,10000);
00305     return true;
00306 }
00307 
00308 void HessianDetector::_cutPointList(double threshold, int nrPoints) {
00309     if(this->orderedList.size()<nrPoints) return;
00310     vector<vector<int> >::iterator iter1 = this->orderedList.begin();
00311     double avg=0;
00312     int count=0;
00313     while( iter1 != this->orderedList.end()) { //loop over every interest point
00314         vector<int > interestPoint=*iter1;
00315 
00316         iter1++;
00317         vector<int > interestPoint2=*iter1;
00318 
00319         //cout << this->determinants[interestPoint[0]][interestPoint[1]] << endl;
00320         //cout << count << endl;
00321 
00322 
00323         iter1--;
00324 
00325         if(this->determinants[interestPoint[0]][interestPoint[1]]<threshold) {
00326                     this->orderedList.erase(iter1);
00327         } else {
00328                 avg+=this->determinants[interestPoint[0]][interestPoint[1]];
00329                 iter1++;
00330         }
00331         count++;
00332     }
00333 
00334     cout << "Number of points:"<< this->orderedList.size()<< endl;
00335 
00336     this->_cutPointList(avg/this->orderedList.size(), nrPoints);
00337 }
00338 
00339 void HessianDetector::_calculateMaxDet(int x, int y) {
00343     int octaves[][4] ={
00344     9,15,21,27,
00345     15,27,39,51,
00346     21,45,69,93
00347     };
00348     double interp1, interp2;
00349 
00350     //first octave
00351     interp1=(this->_convolutePixel(&x,&y,&octaves[0][0])+this->_convolutePixel(&x,&y,&octaves[0][1]))/2;
00352     interp2=(this->_convolutePixel(&x,&y,&octaves[0][2])+this->_convolutePixel(&x,&y,&octaves[0][3]))/2;
00353     if(interp1>interp2) {
00354         this->determinants[x][y]=interp1;
00355         this->maximas[x][y]=(((octaves[0][0]+octaves[0][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
00356     } else {
00357         this->determinants[x][y]=interp2;
00358         this->maximas[x][y]=(((octaves[0][2]+octaves[0][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
00359     }
00360 
00361     //second octave
00362     if(this->nrOctaves>1 && x%2==0 && y%2==0) { //points are sampled
00363         interp1=(this->_convolutePixel(&x,&y,&octaves[1][0])+this->_convolutePixel(&x,&y,&octaves[1][1]))/2;
00364         interp2=(this->_convolutePixel(&x,&y,&octaves[1][2])+this->_convolutePixel(&x,&y,&octaves[1][3]))/2;
00365         if(interp1>this->determinants[x][y] || interp2>this->determinants[x][y] ) {
00366             if(interp1>interp2) {
00367                 this->determinants[x][y]=interp1;
00368                 this->maximas[x][y]=(((octaves[1][0]+octaves[1][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
00369             } else {
00370                 this->determinants[x][y]=interp2;
00371                 this->maximas[x][y]=(((octaves[1][2]+octaves[1][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
00372             }
00373         }
00374     }
00375 
00376     //third octave
00377     if(this->nrOctaves==HD_MAX_OCTAVES && x%4==0 && y%4==0) { //points are sampled
00378         interp1=(this->_convolutePixel(&x,&y,&octaves[2][0])+this->_convolutePixel(&x,&y,&octaves[2][1]))/2;
00379         interp2=(this->_convolutePixel(&x,&y,&octaves[2][2])+this->_convolutePixel(&x,&y,&octaves[2][3]))/2;
00380         if(interp1>this->determinants[x][y] || interp2>this->determinants[x][y] ) {
00381             if(interp1>interp2) {
00382                 this->determinants[x][y]=interp1;
00383                 this->maximas[x][y]=(((octaves[2][0]+octaves[2][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
00384             } else {
00385                 this->determinants[x][y]=interp2;
00386                 this->maximas[x][y]=(((octaves[2][2]+octaves[2][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
00387             }
00388         }
00389     }
00390 }
00391 
00392 int HessianDetector::_getHessianDeterminant(int* pixelSumXX, int* pixelSumXY, int* pixelSumYY) {
00393  return *pixelSumXX*(*pixelSumYY)-pow((0.9*(*pixelSumXY)),2);
00394 }
00395 
00396 int HessianDetector::_convolutePixel(int* x, int* y, int* kernelSize) {
00397 //for box filters only
00398     int pixelSumXX=0;
00399     int pixelSumYY=0;
00400     int pixelSumXY=0;
00401     int det;
00402     switch(*kernelSize) {
00403         case 9:
00404             //9x9 kernel
00405             //getRegionSum(y1,x1,y2,x2)
00406             pixelSumYY=this->image->getRegionSum(*y-4,*x-2,*y-2,*x+2);
00407             pixelSumYY+=this->image->getRegionSum(*y+2,*x-2,*y+4,*x+2);
00408             pixelSumYY+=-2*this->image->getRegionSum(*y-1,*x-2,*y+1,*x+2);
00409 
00410             pixelSumXY=-1*this->image->getRegionSum(*y-3,*x+1,*y-1,*x+3);
00411             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-3,*y+3,*x-1);
00412             pixelSumXY+=this->image->getRegionSum(*y-3,*x-3,*y-1,*x-1);
00413             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+3,*x+3);
00414 
00415             pixelSumXX=this->image->getRegionSum(*y-2,*x-4,*y+2,*x-2);
00416             pixelSumXX+=this->image->getRegionSum(*y-2,*x+2,*y+2,*x+4);
00417             pixelSumXX+=-2*this->image->getRegionSum(*y-2,*x-1,*y+2,*x+1);
00418 
00419             //determinant of the Hessian Matrix
00420                         det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY);
00421         break;
00422         case 11:
00423             //11x11 kernel
00424             pixelSumYY=this->image->getRegionSum(*y-5,*x-3,*y-2,*x+3);
00425             pixelSumYY+=this->image->getRegionSum(*y+2,*x-3,*y+5,*x+3);
00426             pixelSumYY+=-2*this->image->getRegionSum(*y-1,*x-3,*y+1,*x+3);
00427 
00428             pixelSumXY=-1*this->image->getRegionSum(*y-4,*x+1,*y-1,*x+4);
00429             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-4,*y+4,*x-1);
00430             pixelSumXY+=this->image->getRegionSum(*y-4,*x-4,*y-1,*x-1);
00431             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+4,*x+4);
00432 
00433             pixelSumXX=this->image->getRegionSum(*y-3,*x-5,*y+3,*x-2);
00434             pixelSumXX+=this->image->getRegionSum(*y-3,*x+2,*y+3,*x+5);
00435             pixelSumXX+=-2*this->image->getRegionSum(*y-3,*x-1,*y+3,*x+1);
00436             //determinant of the Hessian Matrix
00437                         det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/2.232;
00438             break;
00439         case 15:
00440             //15x15 kernel
00441             pixelSumYY=this->image->getRegionSum(*y-7,*x-4,*y-3,*x+4);
00442             pixelSumYY+=this->image->getRegionSum(*y+3,*x-4,*y+7,*x+4);
00443             pixelSumYY+=-2*this->image->getRegionSum(*y-2,*x-4,*y+2,*x+4);
00444 
00445             pixelSumXY=-1*this->image->getRegionSum(*y-5,*x+1,*y-1,*x+5);
00446             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-5,*y+5,*x-1);
00447             pixelSumXY+=this->image->getRegionSum(*y-5,*x-5,*y-1,*x-1);
00448             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+5,*x+5);
00449 
00450             pixelSumXX=this->image->getRegionSum(*y-4,*x-7,*y+4,*x-3);
00451             pixelSumXX+=this->image->getRegionSum(*y-4,*x+3,*y+4,*x+7);
00452             pixelSumXX+=-2*this->image->getRegionSum(*y-4,*x-2,*y+4,*x+2);
00453 
00454             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/7.72;
00455             break;
00456         case 17:
00457             //17x17 kernel
00458             pixelSumYY=this->image->getRegionSum(*y-8,*x-4,*y-3,*x+4);
00459             pixelSumYY+=this->image->getRegionSum(*y+3,*x-4,*y+8,*x+4);
00460             pixelSumYY+=-2*this->image->getRegionSum(*y-2,*x-4,*y+2,*x+4);
00461 
00462             pixelSumXY=-1*this->image->getRegionSum(*y-6,*x+1,*y-1,*x+6);
00463             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-6,*y+6,*x-1);
00464             pixelSumXY+=this->image->getRegionSum(*y-6,*x-6,*y-1,*x-1);
00465             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+6,*x+6);
00466 
00467             pixelSumXX=this->image->getRegionSum(*y-4,*x-8,*y+4,*x-3);
00468             pixelSumXX+=this->image->getRegionSum(*y-4,*x+3,*y+4,*x+8);
00469             pixelSumXX+=-2*this->image->getRegionSum(*y-4,*x-2,*y+4,*x+2);
00470             //determinant of the Hessian Matrix
00471                         det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/12.73;
00472             break;
00473         case 21:
00474             //21x21 kernel
00475             pixelSumYY=this->image->getRegionSum(*y-10,*x-6,*y-4,*x+6);
00476             pixelSumYY+=this->image->getRegionSum(*y+4,*x-6,*y+10,*x+6);
00477             pixelSumYY+=-2*this->image->getRegionSum(*y-3,*x-6,*y+3,*x+6);
00478 
00479             pixelSumXY=-1*this->image->getRegionSum(*y-7,*x+1,*y-1,*x+7);
00480             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-7,*y+7,*x-1);
00481             pixelSumXY+=this->image->getRegionSum(*y-7,*x-7,*y-1,*x-1);
00482             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+7,*x+7);
00483 
00484             pixelSumXX=this->image->getRegionSum(*y-6,*x-10,*y+6,*x-4);
00485             pixelSumXX+=this->image->getRegionSum(*y-6,*x+4,*y+6,*x+10);
00486             pixelSumXX+=-2*this->image->getRegionSum(*y-6,*x-3,*y+6,*x+3);
00487 
00488             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/29.642;
00489             break;
00490         case 27:
00491             //27x27 kernel
00492             pixelSumYY=this->image->getRegionSum(*y-13,*x-7,*y-5,*x+7);
00493             pixelSumYY+=this->image->getRegionSum(*y+5,*x-7,*y+13,*x+7);
00494             pixelSumYY+=-2*this->image->getRegionSum(*y-4,*x-7,*y+4,*x+7);
00495 
00496             pixelSumXY=-1*this->image->getRegionSum(*y-9,*x+1,*y-1,*x+9);
00497             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-9,*y+9,*x-1);
00498             pixelSumXY+=this->image->getRegionSum(*y-9,*x-9,*y-1,*x-1);
00499             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+9,*x+9);
00500 
00501             pixelSumXX=this->image->getRegionSum(*y-7,*x-13,*y+7,*x-5);
00502             pixelSumXX+=this->image->getRegionSum(*y-7,*x+5,*y+7,*x+13);
00503             pixelSumXX+=-2*this->image->getRegionSum(*y-7,*x-4,*y+7,*x+4);
00504 
00505             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/81;
00506             break;
00507         case 39:
00508             //39x39 kernel
00509             pixelSumYY=this->image->getRegionSum(*y-19,*x-11,*y-7,*x+11);
00510             pixelSumYY+=this->image->getRegionSum(*y+7,*x-11,*y+19,*x+11);
00511             pixelSumYY+=-2*this->image->getRegionSum(*y-6,*x-11,*y+6,*x+11);
00512 
00513             pixelSumXY=-1*this->image->getRegionSum(*y-13,*x+1,*y-1,*x+13);
00514             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-13,*y+13,*x-1);
00515             pixelSumXY+=this->image->getRegionSum(*y-13,*x-13,*y-1,*x-1);
00516             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+13,*x+13);
00517 
00518             pixelSumXX=this->image->getRegionSum(*y-11,*x-19,*y+11,*x-7);
00519             pixelSumXX+=this->image->getRegionSum(*y-11,*x+7,*y+11,*x+19);
00520             pixelSumXX+=-2*this->image->getRegionSum(*y-11,*x-6,*y+11,*x+6);
00521 
00522             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/352.61;
00523             break;
00524         case 45:
00525             //45x45 kernel
00526             pixelSumYY=this->image->getRegionSum(*y-22,*x-12,*y-8,*x+12);
00527             pixelSumYY+=this->image->getRegionSum(*y+8,*x-12,*y+22,*x+12);
00528             pixelSumYY+=-2*this->image->getRegionSum(*y-7,*x-12,*y+7,*x+12);
00529 
00530             pixelSumXY=-1*this->image->getRegionSum(*y-15,*x+1,*y-1,*x+15);
00531             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-15,*y+15,*x-1);
00532             pixelSumXY+=this->image->getRegionSum(*y-15,*x-15,*y-1,*x-1);
00533             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+15,*x+15);
00534 
00535             pixelSumXX=this->image->getRegionSum(*y-12,*x-22,*y+12,*x-8);
00536             pixelSumXX+=this->image->getRegionSum(*y-12,*x+8,*y+12,*x+22);
00537             pixelSumXX+=-2*this->image->getRegionSum(*y-12,*x-7,*y+12,*x+7);
00538 
00539             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/625;
00540             break;
00541         case 51:
00542             //51x51 kernel
00543             pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
00544             pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
00545             pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
00546 
00547             pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
00548             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
00549             pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
00550             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
00551 
00552             pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
00553             pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
00554             pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
00555 
00556             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/1031.123;
00557             break;
00558         case 69:
00559 // TODO (zoran#1#):
00560             //69x69 kernel
00561             pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
00562             pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
00563             pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
00564 
00565             pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
00566             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
00567             pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
00568             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
00569 
00570             pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
00571             pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
00572             pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
00573 
00574             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/32.1;
00575             break;
00576         case 93:
00577 // TODO (zoran#1#):
00578             //93x93 kernel
00579             pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
00580             pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
00581             pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
00582 
00583             pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
00584             pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
00585             pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
00586             pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
00587 
00588             pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
00589             pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
00590             pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
00591 
00592             det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/32.1;
00593             break;
00594     }
00595     return det;
00596 }
00597 
00598 void HessianDetector::_insertToList(int* x, int* y) {
00599     vector<int> point;
00600     point.resize(2);
00601     point[0]=*x;
00602     point[1]=*y;
00603 
00604     if(this->orderedList.size()>=this->nrPoints) {
00605         vector<int> tmp=this->orderedList.back();
00606         if(determinants[*x][*y]<=determinants[tmp[0]][tmp[1]]) return;
00607 
00608         tmp = this->orderedList.front();
00609         if(determinants[*x][*y]>=determinants[tmp[0]][tmp[1]]) {
00610                 this->orderedList.insert(this->orderedList.begin(),point);
00611                 this->orderedList.pop_back();
00612                 return;
00613         }
00614         this->orderedList.pop_back();
00615     } else {
00616         if(this->orderedList.size()==0) {
00617                 this->orderedList.push_back(point);
00618                 return;
00619         }
00620         vector<int> tmp=this->orderedList.back();
00621         if(determinants[*x][*y]<=determinants[tmp[0]][tmp[1]]) {
00622             this->orderedList.push_back(point);
00623             return;
00624         }
00625         tmp = this->orderedList.front();
00626         if(determinants[*x][*y]>=determinants[tmp[0]][tmp[1]]) {
00627                 orderedList.insert(orderedList.begin(),point);
00628                 return;
00629         }
00630     }
00631      vector<vector<int> >::iterator iter1 = orderedList.begin();
00632      while( iter1 != orderedList.end()) {
00633          vector<int > tmp=*iter1;
00634          if(determinants[tmp[0]][tmp[1]]<=determinants[*x][*y]) {
00635                 orderedList.insert(iter1,point);
00636                 return;
00637          }
00638        iter1++;
00639      }
00640 }
00641 void HessianDetector::printPoints() {
00642     vector<vector<int> >::iterator iter1 = orderedList.begin();
00643     //int c=0;
00644      while( iter1 != orderedList.end()) {
00645          vector<int > tmp2=*iter1;
00646             cout << "("<<tmp2[0]<<","<<tmp2[1]<<","<<tmp2[2]<<")\n";
00647                 this->image->drawCircle(tmp2[0],tmp2[1],maximas[tmp2[0]][tmp2[1]]*HD_INIT_KERNEL_SIZE);
00648                 //c++;
00649                 //this->image->drawCircle(tmp2[1],tmp2[0],1);
00650        iter1++;
00651      }
00652       this->image->show();
00653      //cout << c <<"\n";
00654 }
00655 void HessianDetector::printPoints(std::ostream & o) {
00656     cout << "Print HD points"<<endl;
00657 
00658     vector<vector<int> >::iterator iter1 = orderedList.begin();
00659     int n=0;
00660 
00661 /*
00662     while( iter1 != orderedList.end()) {
00663         n++;
00664         iter1++;
00665     }*/
00666     o << 1 << endl;
00667     o << this->orderedList.size() << endl;
00668     iter1 = orderedList.begin();
00669     //int c=0;
00670     while( iter1 != orderedList.end()) {
00671          vector<int > tmp2=*iter1;
00672          //double r = getMaxima(tmp2[0], tmp2[1])*HD_INIT_KERNEL_SIZE;
00673          double r = 1;
00674          o <<tmp2[1]<<" "<<tmp2[0]<<" "<< " " << 1/(r*r) << " " << 0 << " " <<  1/(r*r) << endl;
00675          iter1++;
00676      }
00677     cout << "EOF Print HD points"<<endl;
00678 
00679 }
00680 
00681 vector<vector<int> >* HessianDetector::getPoints() {
00682     return &this->orderedList;
00683 }
00684 int HessianDetector::getNrPoints() {
00685     return this->nrPoints;
00686 }
00687 double HessianDetector::getMaxima(int x, int y) {
00688     return this->maximas[x][y];
00689 }
00690 double HessianDetector::_getScale(int kernelSize) {
00691     return kernelSize/HD_INIT_KERNEL_SIZE;
00692 }
00693 void HessianDetector::dump() {
00694     delete(&determinants);
00695 }

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