Revision: 3154
http://hugin.svn.sourceforge.net/hugin/?rev=3154&view=rev
Author: onurcc
Date: 2008-06-28 04:33:30 -0700 (Sat, 28 Jun 2008)
Log Message:
-----------
Feature matching codes so far. Also added Zoran's feature extractor codes to control_points.
Added Paths:
-----------
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingAlgorithm.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/KDTreeKeypointMatcher.cpp
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/KDTreeKeypointMatcher.h
hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/edgedetection.hxx
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,423 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <string>
+#include <math.h>
+#include "APImage.h"
+
+#ifdef USE_VIGRA
+ using namespace vigra;
+#endif
+
+APImage::APImage(string p) {
+this->path=p;
+}
+
+bool APImage::open() {
+
+#ifdef USE_VIGRA
+ try
+ {
+ // read image given as first argument
+ // file type is determined automatically
+ ImageImportInfo info(this->path.c_str());
+ //cout << this->path.c_str() << "Size:"<< info.width() <<"\n";
+
+ if(info.isGrayscale())
+ {
+ // create a gray scale image of appropriate size
+ vigra::BImage in(info.width(), info.height());
+ this->imgBW=new BImage(info.width(), info.height());
+ // import the image just read
+ importImage(info, destImage(*this->imgBW));
+ }
+ else
+ {
+
+ BRGBImage in(info.width(), info.height());
+
+ //TODO: this does not work. How to make a conversion from RGB to BW image in vigra???
+ this->imgBW=new BImage(info.width(), info.height());
+
+ importImage(info, destImage(in));
+
+ // create image iterator that points to upper left corner
+ // of source image
+ vigra::BRGBImage::Iterator sy = in.upperLeft();
+
+ // create image iterator that points past the lower right corner of
+ // source image (similarly to the past-the-end iterator in the STL)
+ vigra::BRGBImage::Iterator send = in.lowerRight();
+
+ // create image iterator that points to upper left corner
+ // of destination image
+ vigra::BImage::Iterator dy = (*this->imgBW).upperLeft();
+
+ // iterate down the first column of the images
+ for(; sy.y != send.y; ++sy.y, ++dy.y)
+ {
+ // create image iterator that points to the first
+ // pixel of the current row of the source image
+ vigra::BRGBImage::Iterator sx = sy;
+
+ // create image iterator that points to the first
+ // pixel of the current row of the destination image
+ vigra::BImage::Iterator dx = dy;
+
+ // iterate across current row
+ for(; sx.x != send.x; ++sx.x, ++dx.x)
+ {
+ RGBValue<int,0u,1u,2u> pixel = *sx;
+ // calculate grayscale value
+ // Y = 0.3*R + 0.59*G + 0.11*B
+ *dx = vigra::round(0.3*pixel.red() + 0.59*pixel.green() + 0.11*pixel.blue());
+ }
+ }
+ //exportAPImage(srcAPImageRange(*this->imgBW), vigra::APImageExportInfo("ttt.jpg"));
+ }
+ return true;
+ }
+ catch (vigra::StdException & e)
+ {
+ // catch any errors that might have occured and print their reason
+ return false;
+ }
+#endif
+}
+void APImage::convolute(int* kernel,int dim1, int dim2,double scale) {
+ /*CvScalar s;
+ int offsetX=floor(dim1/2);
+ int offsetY=floor(dim2/2);
+ int xStart=0;
+ int yStart=0;
+
+ double pixelSum=0;
+ // bool outOfRange=false;
+ int tmpX=0;
+ int tmpY=0;
+ double max=0;
+
+ for(int i=0;i<imgBW->height;i++) {
+ for(int j=0; j<imgBW->width;j++) {
+ //s=cvGet2D(imgBW,i,j);
+ outOfRange=true;
+ if(j-offsetX<0) {
+ outOfRange=true;
+ }
+ if(i-offsetX<0)
+
+ pixelSum=0;
+ xStart=i-offsetX;
+ yStart=j-offsetY;
+ for(int k=0;k<dim1;k++) {
+ tmpX=xStart+k;
+ for(int l=0;l<dim2;l++) {
+ tmpY=yStart+l;
+ if((tmpX)<0 ||tmpY<0 || tmpX>=imgBW->height || tmpY>=imgBW->width) {
+ pixelSum+=kernel[k*dim2+l]*(rand()%256);
+ } else {
+
+ s=cvGet2D(imgBW,xStart+k,yStart+l);
+ //cout << " "<< s.val[0] << "\n";
+
+ pixelSum+=kernel[k*dim2+l]*s.val[0];
+ }
+ //cout << "\nPixelSum:"<< pixelSum;
+ }
+ }
+ s.val[0]=pixelSum*0.004;
+ cvSet2D(imgBW,i,j,s); // set the (i,j) pixel value
+ if(pixelSum>max) max=pixelSum;
+ //cout << " "<< pixelSum;
+ //this->convolution[i][j]=pixelSum;
+ //s=cvGet2D(imgBW,i,j);
+// s.val[0]=rand()%255;
+// cvSet2D(imgBW,i,j,s); // set the (i,j) pixel value
+ }
+ //out << "\n";
+ }*/
+}
+int APImage::getPixel(int x,int y) {
+#ifdef USE_VIGRA
+ return (*this->imgBW)(x,y);
+#endif
+#ifdef USE_OPENCV
+ CvScalar s=cvGet2D(imgBW,x,y);
+ s.val[0];
+#endif
+}
+int APImage::getWidth() {
+ return this->imgBW->width();
+}
+int APImage::getWidthBW() {
+ return this->imgBW->width();
+}
+int APImage::getHeight() {
+ return this->imgBW->height();
+}
+int APImage::getHeightBW() {
+ return this->imgBW->height();
+}
+APImage* APImage::getCopy() {
+ return new APImage("");
+}
+void APImage::scale(double factor) {
+//#ifdef USE_OPENCV
+ /*IplAPImage *resized= cvCreateAPImage(cvSize(round(this->getWidthBW()/factor),round(this->getHeightBW()/factor)), imgBW->depth ,imgBW->nChannels);
+ cvResize(imgBW,resized);
+ this->imgBW=resized;*/
+//#endif
+}
+void APImage::show() {
+#ifdef USE_VIGRA
+ string add = "00det_";
+ add.append(this->getPath());
+
+ cout << "Results:"<< add << endl;
+ exportImage(srcImageRange(*this->imgBW), ImageExportInfo(add.c_str()));
+#endif
+#ifdef USE_OPENCV
+ cvNamedWindow( "Image view", 1 );
+ cvShowAPImage( "Image view", this->img);
+ cvWaitKey(0); // very important, contains event processing loop inside
+ cvDestroyWindow( "Image view" );
+#endif
+ //cvReleaseAPImage( &img );
+}
+void APImage::drawCircle(int x,int y, int radius) {
+ //cout << "Circle:" << x << "," << y << ",radius:"<<radius<<"\n";
+#ifdef USE_VIGRA
+
+ (*this->imgBW)(x,y) = 255;
+ (*this->imgBW)(x+1,y+1) = 255;
+ (*this->imgBW)(x-1,y-1) = 255;
+ (*this->imgBW)(x-1,y+1) = 255;
+ (*this->imgBW)(x+1,y-1) = 255;
+
+#endif
+ //cout << "EOF CIRCLE DRAWING" <<endl;
+
+#ifdef USE_OPENCV
+ cvCircle(this->img, cvPoint(x,y), radius, cvScalar(0,255,0), 1);
+ cvCircle(this->img, cvPoint(x,y), 0, cvScalar(0,255,0), 1);
+#endif
+
+}
+void APImage::drawRectangle(int x,int y, int radius) {
+ //cvRectangle(this->img,cvPoint(x-radius,y-radius),cvPoint(x+radius,y+radius),cvScalar(0,255,0));
+}
+
+void APImage::drawLine(int x1,int y1, int x2,int y2) {
+ //cout << "Draw line:("<< x1<< ","<<y1<<")->("<<x2<<","<<y2<<")\n";
+#ifdef USE_OPENCV
+ cvLine(this->img,cvPoint(x1,y1),cvPoint(x2,y2),cvScalar(0,255,255));
+#endif
+
+}
+
+void APImage::smooth() {
+ //TODO: create a gaussian mask and call the function convolute
+}
+/**
+ * Calculates the integral image
+ */
+void APImage::integrate() {
+
+ /*cout << "Height:"<< this->getHeightBW() <<"\n";
+ cout << "Width:"<< this->getWidthBW() <<"\n";*/
+
+ /*for(int i=0;i<10;i++) {
+ for(int j=0; j<10;j++) {
+ cout << this->getPixel(j,i)<<" ";
+ }
+ cout << "\n";
+ }*/
+
+ this->integral.clear();
+ this->integral.resize(this->getWidthBW());
+ for(int i=0;i<this->getWidthBW();i++) {
+ this->integral[i].resize(this->getHeightBW());
+ for(int j=0; j<this->getHeightBW();j++) {
+ //cout << i << ","<< j<<"\n";
+ this->integral[i][j]=this->_getValue4Integral(i,j-1)+this->_getValue4Integral(i-1,j)+this->getPixel(i,j)-this->_getValue4Integral(i-1,j-1);
+ }
+ //cout << "\n";
+ }
+
+ /*for(int i=0;i<10;i++) {
+ for(int j=0; j<10;j++) {
+ cout << this->getIntegralPixel(j,i)<<" ";
+ }
+ cout << "\n";
+ }
+
+ cout << this->getRegionSum(2,1,2,3) <<"dddddd\n";*/
+}
+int APImage::_getValue4Integral(int x, int y) {
+if(x==-1 || y==-1) return 0;
+else return this->integral[x][y];
+}
+int APImage::getIntegralPixel(int x,int y) {
+ return this->integral[x][y];
+}
+int APImage::getRegionSum(int y1, int x1, int y2, int x2) {
+ if(x1<=0) x1=1;
+ if(y1<=0) y1=1;
+ if(x2<=0) x2=1;
+ if(y2<=0) y2=1;
+ if(x1>=this->getWidthBW()) x1=this->getWidthBW()-1;
+ if(x2>=this->getWidthBW()) x2=this->getWidthBW()-1;
+ if(y1>=this->getHeightBW()) y1=this->getHeightBW()-1;
+ if(y2>=this->getHeightBW()) y2=this->getHeightBW()-1;
+ return this->integral[x2][y2]+this->integral[x1-1][y1-1]-this->integral[x1-1][y2]-this->integral[x2][y1-1];
+}
+string APImage::getPath() {
+ return this->path;
+}
+
+/*double APImage::getPointOrientation(int x, int y) {
+
+ return 0.0;
+}*/
+
+template <class SrcIterator, class SrcAccessor, class BackInsertable>
+void APImage::_cannyEdgelList1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
+ BackInsertable & edgels, double scale,vector<int>* point )
+{
+ _cannyEdgelList(src.first, src.second, src.third, edgels, scale, point);
+}
+
+template <class SrcIterator, class SrcAccessor, class BackInsertable>
+void APImage::_cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
+ BackInsertable & edgels, double scale, vector<int>* point)
+{
+ int w = lr.x - ul.x;
+ int h = lr.y - ul.y;
+
+ // calculate image gradients
+ typedef typename
+ NumericTraits<typename SrcAccessor::value_type>::RealPromote
+ TmpType;
+
+ BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h);
+
+ Kernel1D<double> smooth, grad;
+
+ smooth.initGaussian(scale);
+ grad.initGaussianDerivative(scale, 1);
+
+ separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
+ separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth));
+
+ separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
+ separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth));
+
+ combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp),
+ MagnitudeFunctor<TmpType>());
+
+
+ // find edgels
+ internalCannyFindEdgels(dx, dy, tmp, edgels, point);
+}
+
+template <class Image1, class Image2, class BackInsertable>
+void APImage::_internalCannyFindEdgels(Image1 const & gx,
+ Image1 const & gy,
+ Image2 const & magnitude,
+ BackInsertable & edgels, vector<int>* p)
+{
+ typedef typename Image1::value_type PixelType;
+ double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
+
+ //last element in edgel list is edgel that holds orientation
+ //of interest point
+
+ //orientation assignment
+ vector<int > point= *p;
+ PixelType gradx = gx(p[0],p[1]);
+ PixelType grady = gy(p[0],p[1]);
+
+ double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
+ if(orientation < 0.0)
+ orientation += 2.0*M_PI;
+ Edgel edgel;
+ edgel.orientation=orientation;
+ edgels.push_back(edgel);
+ //EOF orientation assignment
+
+ for(int y=1; y<gx.height()-1; ++y)
+ {
+ for(int x=1; x<gx.width()-1; ++x)
+ {
+ gradx = gx(x,y);
+ grady = gy(x,y);
+ double mag = magnitude(x, y);
+
+ int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
+ int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
+
+ int x1 = x - dx,
+ x2 = x + dx,
+ y1 = y - dy,
+ y2 = y + dy;
+
+ PixelType m1 = magnitude(x1, y1);
+ PixelType m3 = magnitude(x2, y2);
+
+ if(m1 < mag && m3 <= mag)
+ {
+ Edgel edgel;
+
+ // local maximum => quadratic interpolation of sub-pixel location
+ PixelType del = (m1 - m3) / 2.0 / (m1 + m3 - 2.0*mag);
+ edgel.x = x + dx*del;
+ edgel.y = y + dy*del;
+ edgel.strength = mag;
+ orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
+ if(orientation < 0.0)
+ orientation += 2.0*M_PI;
+ edgel.orientation = orientation;
+ edgels.push_back(edgel);
+ }
+ }
+ }
+}
+
+void APImage::test() {
+ // empty edgel list
+ /*std::vector<vigra::Edgel> edgels;
+
+ // find edgels at scale of the interest point
+ vigra::cannyEdgelList(srcIterRange( this->imgBW->upperLeft() + vigra::Diff2D(250, 250),
+ this->imgBW->upperLeft() + vigra::Diff2D(550, 450)),
+ edgels, 1.2, );
+ cout << "Size:" << edgels.size() << endl;
+ vector<vigra::Edgel>::iterator iter2 = edgels.begin();
+ int a=0;
+ while( iter2 != edgels.end()) { //loop over every canny pixel
+ vigra::Edgel edgePoint=*iter2;
+ //cout << edgePoint.strength << endl;
+ if(edgePoint.strength>5) {
+ this->drawCircle(round(edgePoint.x)+250,round(edgePoint.y)+250,0);
+ a++;
+ }
+ iter2++;
+ }
+ cout << a << endl;*/
+}
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/APImage.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,121 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#ifndef IMAGE_H_INCLUDED
+#define IMAGE_H_INCLUDED
+#define USE_VIGRA
+
+#include <string>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <vector>
+
+#ifdef USE_OPENCV
+#include "cv.h"
+#endif
+
+#ifdef USE_VIGRA
+#include "vigra/stdimage.hxx"
+#include "vigra/imageinfo.hxx"
+#include "vigra/impex.hxx"
+#include "vigra/stdimagefunctions.hxx"
+#include "edgedetection.hxx"
+#include "vigra/utilities.hxx"
+#include "vigra/numerictraits.hxx"
+
+#include "vigra/recursiveconvolution.hxx"
+#include "vigra/separableconvolution.hxx"
+#include "vigra/labelimage.hxx"
+#include "vigra/mathutil.hxx"
+#include "vigra/pixelneighborhood.hxx"
+#include "vigra/linear_solve.hxx"
+#endif
+
+
+using namespace std;
+ using namespace vigra;
+
+class APImage
+ {
+
+ public:
+ APImage(string p);
+ void convolute(int* kernel,int dim1, int dim2,double scale);
+ string getPath();
+ int getWidth();
+ int getWidthBW();
+ int getHeight();
+ int getHeightBW();
+ int getPixel(int x, int y);
+ int getIntegralPixel(int x,int y);
+ void scale(double factor);
+ APImage* getCopy();
+ void drawCircle(int x,int y, int radius);
+ void drawLine(int x1,int y1, int x2,int y2);
+ void drawRectangle(int x,int y, int radius);
+ void smooth();
+ void integrate();
+ int getRegionSum(int x1, int y1, int x2, int y2);
+
+ bool open();
+ void show();
+ void test();
+
+ template <class SrcIterator, class SrcAccessor, class BackInsertable>
+ void _cannyEdgelList1(vigra::triple<SrcIterator, SrcIterator, SrcAccessor> src,
+ BackInsertable & edgels, double scale, vector<int>* p);
+
+
+#ifdef USE_VIGRA
+ vigra::BImage* imgBW;
+#endif
+
+/* private slots:
+*/
+
+ private:
+ /**
+ * Holds the convolution of the image.
+ */
+ vector<vector<int> > convolution;
+ /**
+ * Holds the values of the integral image.
+ */
+ vector<vector<int> > integral;
+ string path;
+ int _getValue4Integral(int x, int y);
+
+ template <class SrcIterator, class SrcAccessor, class BackInsertable>
+ inline void _cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
+ BackInsertable & edgels, double scale, vector<int>* point);
+
+ template <class Image1, class Image2, class BackInsertable>
+ void _internalCannyFindEdgels(Image1 const & gx,
+ Image1 const & gy,
+ Image2 const & magnitude,
+ BackInsertable & edgels, vector<int>* p);
+
+ /*IplImage* img;
+ IplImage* imgBW;*/
+
+ };
+
+
+#endif // IMAGE_H_INCLUDED
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,839 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#include <stdio.h>
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+#include "Descriptor.h"
+
+using namespace std;
+
+Descriptor::Descriptor(APImage* i,HessianDetector* hessianDetector) {
+ this->image=i;
+ this->hd=hessianDetector;
+}
+
+void Descriptor::setPoints(vector<vector<int> >* pts) {
+ this->interestPoints=pts;
+}
+
+// dangelo: this does not compile with MSVC and seems to be unused anyway...
+#if 0
+void Descriptor::orientate() {
+ vector<vector<int> >::iterator iter1 = this->interestPoints->begin();
+ double regionSize;
+ int kernelSize;
+
+ //int angle;
+ int max;
+ int orientation;
+ double step;
+
+ //coordinates of the sample
+ int pointX;
+ int pointY;
+
+/*
+S0: if dx>0 and dy>0 and |dy|>|dx|
+S1: if dx>0 and dy>0 and |dy|<=|dx|
+S2: if dx>0 and dy<=0 and |dy|<=|dx|
+S3: if dx>0 and dy<=0 and |dy|>|dx|
+S4: if dx<=0 and dy<=0 and |dy|>|dx|
+S5: if dx<=0 and dy<=0 and |dy|<=|dx|
+S6: if dx<=0 and dy>0 and |dy|<=|dx|
+S7: if dx<=0 and dy>0 and |dy|>|dx|
+ */
+ double Sx[NR_ANGLE_BINS];
+ double Sy[NR_ANGLE_BINS];
+
+
+
+ double dx;
+ double dy;
+
+ int pointCount=0;
+
+ while( iter1 != this->interestPoints->end()) { //loop over every interest point
+ vector<int > interestPoint=*iter1;
+ //cout << interestPoint[0]<<","<<interestPoint[1]<<":";
+ regionSize=_getMaxima(interestPoint[0], interestPoint[1])*3;
+ kernelSize=vigra::round(_getMaxima(interestPoint[0], interestPoint[1])*2);
+ step=_getMaxima(interestPoint[0], interestPoint[1]);
+
+ for(int i=0; i<NR_ANGLE_BINS;i++) {
+ Sx[i]=0;
+ Sy[i]=0;
+ }
+
+ /*for(int i=interestPoint[0]-regionSize;i<=interestPoint[0]+regionSize;i++) {
+ for(int j=interestPoint[1]-regionSize; j<=interestPoint[1]+regionSize;j++) {
+ cout << "(" <<this->image->getPixel(i,j)<< "," << this->image->getIntegralPixel(i,j) <<") ";
+ }
+ cout << "\n";
+ }*/
+
+ for(double i=(interestPoint[0]-regionSize);i<=(interestPoint[0]+regionSize);i+=step) {
+ for(double j=(interestPoint[1]-regionSize); j<=(interestPoint[1]+regionSize);j+=step) {
+
+ //pixels have integer values
+ pointX=vigra::round(i);
+ pointY=vigra::round(j);
+ //pointX and pointY are the coordinates of the sample point
+
+ double dist=_euclidianDistance(interestPoint[0],interestPoint[1], pointX,pointY);
+ if(dist<=regionSize) { //if point is in the circle
+
+ double weight=_gaussWeighting(interestPoint[0]-pointX,interestPoint[1]-pointY,2.5*_getMaxima(pointX,pointY));
+
+ //x directon
+ dx=this->image->getRegionSum(pointX-kernelSize, pointY,pointX+kernelSize,pointY+kernelSize);
+ dx-=this->image->getRegionSum(pointX-kernelSize, pointY-kernelSize,pointX+kernelSize,pointY);
+ dx*=weight;
+
+ //y direction
+ dy=this->image->getRegionSum(pointX,pointY-kernelSize,pointX+kernelSize, pointY+kernelSize);
+ dy-=this->image->getRegionSum(pointX-kernelSize,pointY-kernelSize,pointX, pointY+kernelSize);
+ dy*=weight;
+
+ if(dx>0) {
+ if(dy>0) {
+ if(fabs(dy)>fabs(dx)) {
+ Sx[0]+=dx;
+ Sy[0]+=dy;
+ } else {
+ Sx[1]+=dx;
+ Sy[1]+=dy;
+ }
+ } else {
+ if(fabs(dy)<=fabs(dx)) {
+ Sx[2]+=dx;
+ Sy[2]+=dy;
+ } else {
+ Sx[3]+=dx;
+ Sy[3]+=dy;
+ }
+ }
+ } else {
+ if(dy<=0) {
+ if(fabs(dy)>fabs(dx)) {
+ Sx[4]+=dx;
+ Sy[4]+=dy;
+ } else {
+ Sx[5]+=dx;
+ Sy[5]+=dy;
+ }
+ } else {
+ if(fabs(dy)<=fabs(dx)) {
+ Sx[6]+=dx;
+ Sy[6]+=dy;
+ } else {
+ Sx[7]+=dx;
+ Sy[7]+=dy;
+ }
+ }
+ }
+
+ } //if point in circle
+ } //for j
+ } //for i
+ max=0;
+ orientation=0;
+
+ for(int i=0; i<NR_ANGLE_BINS;i++) {
+ double tmp = _euclidianDistance(0,0,Sx[i],Sy[i]);
+ if(tmp>max) {
+ max = tmp;
+ orientation=i;
+ }
+ }
+
+ orientation=vigra::round(atan2(Sy[orientation],Sx[orientation]) * 180 / PI);
+ //if(orientation<0) orientation=360+orientation;
+
+ this->orientations[pointCount]=orientation;
+
+ this->image->drawCircle(interestPoint[1],interestPoint[0],1);
+ this->image->drawLine(interestPoint[1],interestPoint[0],interestPoint[1]+vigra::round(sin((double)orientation)*10),interestPoint[0]+vigra::round(cos((double)orientation)*10));
+ //cout << orientation << "\n";
+ iter1++;
+ //pointCount++;
+ }
+ this->image->show();
+}
+
+#endif
+
+double Descriptor::_gaussWeighting(int x, int y, double stdev) {
+return (1/(2*PI*pow(stdev,2.0)))*exp(-(pow((double)x,2.0)+pow((double)y,2.0))/(2*stdev));
+}
+double Descriptor::_getMaxima(int x,int y) {
+ return this->hd->getMaxima(x,y);
+ //return 1.6;
+}
+double Descriptor::_euclidianDistance(int x1, int y1, int x2, int y2) {
+ return sqrt(pow((double)(x1-x2),2.0)+pow((double)(y1-y2),2.0));
+}
+void Descriptor::_GaborResponse(int x,int y, double maxima, double* descriptor) {
+ //
+ double factor = maxima/HD_INITIAL_SCALE;
+ /*vector<int> d=*descriptor;
+ d[0]++;*/
+
+/**
+ * Gabor filter bank consists of 8 kernels. There are 4 orientations and 2 frequencies.
+ * Orientations: 0, 45, 90 and 135 degrees
+ * First frequency: bandwith 2, phase PI/2, aspect ratio 0.1, wavelength 5
+ * Second frequency: bandwith 4, phase PI/2, aspect ratio 0.1, wavelength 10
+ * Third frequency: bandwith 6, phase PI/2, aspect ratio 0.1, wavelength 8
+ */
+ //relative coordinates(from the center of the kernel) and weights
+ //for points of local maxima for first gabor filter(x,y,weight)
+ double gabor11[6][3]={
+ {-1,0,0.839304},
+ {1,0,-0.839304},
+ {-3,0,-0.190826},
+ {3,0,0.190826},
+ {-6,0,0.0105653},
+ {6,0,-0.0105653}
+ };
+ //second...orientation=0, wavelength=10
+ double gabor12[6][3]={
+ {-2,0,0.839304},
+ {2,0,-0.839304},
+ {-7,0,-0.20568},
+ {7,0,0.20568},
+ {-11,0,0.0133981},
+ {11,0,-0.0133981}
+ };
+
+ double gabor13[10][3]={
+ {-2,0,0.945959},
+ {2,0,-0.945959},
+ {-6,0,-0.606531},
+ {6,0,0.606531},
+ {-10,0,0.249352},
+ {10,0,-0.249352},
+ {-13,0,0.0676238},
+ {13,0,-0.0676238},
+ {-17,0,0.0127725},
+ {17,0,-0.0127725}
+ };
+
+ //orientation=PI/4, wavelength=5
+ double gabor21[6][3]={
+ {-1,-1,0.762278},
+ {-3,-2,-0.201919},
+ {2,3,0.201919},
+ {1,1,0.762278},
+ {3,5,0.0134254},
+ {-5,-3,-0.0134254}
+ };
+
+ //orientation=PI/4, wavelength=10
+ double gabor22[6][3]={
+ {-2,-1,0.844207},
+ {2,1,-0.844207},
+ {-5,-4,-0.213173},
+ {5,4,0.213173},
+ {-9,-7,0.013459},
+ {9,7,-0.013459}
+ };
+
+ double gabor23[10][3]={
+ {-2,-1,0.935087},
+ {2,1,-0.935087},
+ {-4,-4,-0.618035},
+ {4,4,0.618035},
+ {-7,-7,0.255577},
+ {7,7,-0.255577},
+ {-9,-10,-0.0736175},
+ {10,9,0.0736175},
+ {-12,-12,0.0126482},
+ {12,12,-0.0126482}
+ };
+
+ double gabor31[6][3]={
+ {0,-1,0.839304},
+ {0,1,-0.839304},
+ {0,-3,-0.190826},
+ {0,3,0.190826},
+ {0,6,-0.0105653},
+ {0,-6,0.0105653}
+ };
+ double gabor32[6][3]={
+ {0,2,-0.839304},
+ {0,7,0.20568},
+ {0,11,-0.0133981},
+ {0,-2,0.839304},
+ {0,-7,-0.20568},
+ {0,-11,0.0133981}
+ };
+ double gabor33[10][3]={
+ {0,-2,0.945959},
+ {0,2,-0.945959},
+ {0,-6,-0.606531},
+ {0,6,0.606531},
+ {0,-10,0.249352},
+ {0,10,-0.249352},
+ {0,-13,-0.0676238},
+ {0,13,0.0676238},
+ {0,-17,-0.0127725},
+ {0,17,0.0127725}
+ };
+
+ double gabor41[6][3]={
+ {1,-1,0.762278},
+ {-2,+3,0.201919},
+ {2,-3,-0.201919},
+ {-1,1,-0.762278},
+ {3,-5,0.0134254},
+ {-3,5,-0.0134254}
+ };
+ double gabor42[6][3]={
+ {-1,2,-0.844207},
+ {2,-1,0.844207},
+ {4,-5,-0.213173},
+ {-5,4,0.213173},
+ {8,-8,0.013459},
+ {-8,8,-0.013459}
+ };
+ double gabor43[10][3]={
+ {2,-1,0.935087},
+ {-1,2,-0.935087},
+ {4,-4,-0.618035},
+ {-4,4,0.618035},
+ {7,-7,0.255577},
+ {-7,7,-0.255577},
+ {9,-10,-0.0736175},
+ {-10,9,0.0736175},
+ {12,-12,0.0126482},
+ {-12,12,-0.0126482}
+ };
+
+ //orientation PI/8
+ double gabor53[10][3]={
+ {-2,0,0.946801},
+ {2,0,-0.946801},
+ {-5,-3,-0.619484},
+ {3,5,0.619484},
+ {-7,-8,0.263348},
+ {8,7,-0.263348},
+ {-12,-6,0.0735327},
+ {6,12,-0.0735327},
+ {-15,-9,0.0133379},
+ {9,15,-0.0133379}
+ };
+
+ //orientation 3*(PI/8)
+ double gabor63[10][3]={
+ {0,-2,0.946801},
+ {0,2,-0.946801},
+ {-3,-5,-0.619484},
+ {5,3,0.619484},
+ {-8,-7,0.263348},
+ {7,8,-0.263348},
+ {-6,-12,0.0735327},
+ {12,6,-0.0735327},
+ {-9,-15,0.0133379},
+ {15,9,-0.0133379}
+ };
+
+ //orientation 5*PI/8
+ double gabor73[10][3]={
+ {0,-2,0.946801},
+ {0,2,-0.946801},
+ {3,-5,-0.619484},
+ {-3,5,0.619484},
+ {8,-7,0.263348},
+ {-8,8,-0.263348},
+ {6,-12,0.0735327},
+ {-6,12,-0.0735327},
+ {9,-15,0.0133379},
+ {-9,15,-0.0133379}
+ };
+
+ //orientation 7*PI/8
+ double gabor83[10][3]={
+ {2,0,0.946801},
+ {-2,0,-0.946801},
+ {5,-3,-0.619484},
+ {-5,3,0.619484},
+ {7,-8,0.263348},
+ {-7,8,-0.263348},
+ {12,-6,0.0735327},
+ {-12,6,-0.0735327},
+ {15,-9,0.0133379},
+ {-15,9,-0.0133379}
+ };
+
+ //stdev=6, wavelength=4, phase=0, orientation=0
+ double gabor91[17][3]={
+ {-16,0, 0.0285655},
+ {-14,0, -0.0657285},
+ {-12,0, 0.135335},
+ {-10,0, -0.249352},
+ {-8,0, 0.411112},
+ {-6,0, -0.606531},
+ {-4,0, 0.800737},
+ {-2,0, -0.945959},
+ {0,0, 1},
+ {16,0, 0.0285655},
+ {14,0, -0.0657285},
+ {12,0, 0.135335},
+ {10,0, -0.249352},
+ {8,0, 0.411112},
+ {6,0, -0.606531},
+ {4,0, 0.800737},
+ {2,0, -0.945959}
+ };
+
+ //stdev=6, wavelength=4, phase=0, orientation=PI/8
+ double gabor101[17][3]={
+ {-16,-3,0.0292416},
+ {-14,-2, -0.0656082},
+ {-12,-2, 0.138167},
+ {-10,-2, -0.248923},
+ {-8,-1, 0.404746},
+ {-6,-1, -0.609707},
+ {-4,-1, 0.787721},
+ {-2,0, -0.926472},
+ {0,0, 1},
+ {16,3,0.0292416},
+ {14,2, -0.0656082},
+ {12,2, 0.138167},
+ {10,2, -0.248923},
+ {8,1, 0.404746},
+ {6,1, -0.609707},
+ {4,1, 0.787721},
+ {2,0, -0.926472}
+ };
+ //stdev=6, wavelength=4, phase=0, orientation=PI/4
+ double gabor111[17][3]={
+ {-11,-11, 0.026607},
+ {-10,-10, -0.0606333},
+ {-9,-8, 0.134318},
+ {-7,-7, -0.253187},
+ {-6,-5, 0.405626},
+ {-4,-4, -0.550271},
+ {-3,-3, 0.722915},
+ {-2,-1, -0.922342},
+ {0,0, 1},
+ {11,11, 0.026607},
+ {10,10, -0.0606333},
+ {8,9, 0.134318},
+ {7,7, -0.253187},
+ {5,6, 0.405626},
+ {4,4, -0.550271},
+ {3,3, 0.722915},
+ {1,2, -0.922342}
+ };
+
+ //3PI/8
+ double gabor121[17][3]={
+ {-5,-15, 0.0295778},
+ {-5,-13, -0.0672137},
+ {-2,-12, 0.138167},
+ {-4,-9, -0.252577},
+ {-1,-8, 0.404746},
+ {-3,-5, -0.588397},
+ {-1,-4, 0.787721},
+ {-0,-2, -0.926472},
+ {0,0, 1},
+ {5,15, 0.0295778},
+ {5,13, -0.0672137},
+ {2,12, 0.138167},
+ {4,9, -0.252577},
+ {1,8, 0.404746},
+ {3,5, -0.588397},
+ {1,4, 0.787721},
+ {0,2, -0.926472}
+ };
+
+ double gabor131[17][3]={
+ {0,-16, 0.0285655},
+ {0,-14, -0.0657285},
+ {0,-12, 0.135335},
+ {0,-10, -0.249352},
+ {0,-8, 0.411112},
+ {0,-6, -0.606531},
+ {0,-4, 0.800737},
+ {0,-2, -0.945959},
+ {0,0, 1},
+ {0,16, 0.0285655},
+ {0,14, -0.0657285},
+ {0,12, 0.135335},
+ {0,10, -0.249352},
+ {0,8, 0.411112},
+ {0,6, -0.606531},
+ {0,4, 0.800737},
+ {0,2, -0.945959}
+ };
+
+ //5PI/8
+ double gabor141[17][3]={
+ {10,-13, 0.0296244},
+ {5,-14, -0.0672137},
+ {7,-10, 0.137794},
+ {4,-9, -0.252577},
+ {6,-6, 0.411811},
+ {1,-6, -0.606531},
+ {3,-3, 0.801129},
+ {3,-1, -0.93537},
+ {0,0, 1},
+ {-10,13, 0.0296244},
+ {-5,14, -0.0672137},
+ {-7,10, 0.137794},
+ {-4,9, -0.252577},
+ {-6,6, 0.411811},
+ {-1,6, -0.606531},
+ {-3,-3, 0.801129},
+ {-3,1, -0.93537}
+ };
+
+ //6PI/8
+ double gabor151[17][3]={
+ {11,-11, 0.026607},
+ {10,-10, -0.0606333},
+ {9,-8, 0.134318},
+ {7,-7, -0.253187},
+ {6,-5, 0.405626},
+ {4,-4, -0.550271},
+ {3,-3, 0.722915},
+ {2,-1, -0.922342},
+ {0,0, 1},
+ {-11,11, 0.026607},
+ {-10,10, -0.0606333},
+ {-8,9, 0.134318},
+ {-7,7, -0.253187},
+ {-5,6, 0.405626},
+ {-4,4, -0.550271},
+ {-3,3, 0.722915},
+ {-1,2, -0.922342}
+ };
+
+ double gabor161[17][3] = {
+ {15,-5, 0.0295778},
+ {13,-5, -0.0672137},
+ {12,-2, 0.138167},
+ {9,-4, -0.252577},
+ {8,-1, 0.404746},
+ {5,-3, -0.588397},
+ {4,-1, 0.787721},
+ {0,-2, -0.926472},
+ {0,0, 1},
+ {-15,5, 0.0295778},
+ {-13,5, -0.0672137},
+ {-12,2, 0.138167},
+ {-9,4, -0.252577},
+ {-8,1, 0.404746},
+ {-5,3, -0.588397},
+ {-4,1, 0.787721},
+ {-2,0, -0.926472}
+ };
+
+
+ for(int i=0; i<6; i++) {
+
+ descriptor[0]+=this->image->getPixel(vigra::round(x+factor*gabor11[i][0]),vigra::round(y+factor*gabor11[i][1]))*gabor11[i][2];
+ descriptor[2]+=this->image->getPixel(vigra::round(x+factor*gabor21[i][0]),vigra::round(y+factor*gabor21[i][1]))*gabor21[i][2];
+ descriptor[4]+=this->image->getPixel(vigra::round(x+factor*gabor31[i][0]),vigra::round(y+factor*gabor31[i][1]))*gabor31[i][2];
+ descriptor[6]+=this->image->getPixel(vigra::round(x+factor*gabor41[i][0]),vigra::round(y+factor*gabor41[i][1]))*gabor41[i][2];
+
+ descriptor[1]+=this->image->getPixel(vigra::round(x+factor*gabor12[i][0]),vigra::round(x+factor*gabor12[i][1]))*gabor12[i][2];
+ descriptor[3]+=this->image->getPixel(vigra::round(x+factor*gabor22[i][0]),vigra::round(x+factor*gabor22[i][1]))*gabor22[i][2];
+ descriptor[5]+=this->image->getPixel(vigra::round(x+factor*gabor32[i][0]),vigra::round(x+factor*gabor32[i][1]))*gabor32[i][2];
+ descriptor[7]+=this->image->getPixel(vigra::round(x+factor*gabor42[i][0]),vigra::round(x+factor*gabor42[i][1]))*gabor42[i][2];
+
+ }
+ for(int i=0; i<10;i++) {
+
+ descriptor[8]+=this->image->getPixel(vigra::round(x+factor*gabor13[i][0]),vigra::round(y+factor*gabor13[i][1]))*gabor13[i][2];
+ descriptor[9]+=this->image->getPixel(vigra::round(x+factor*gabor23[i][0]),vigra::round(y+factor*gabor23[i][1]))*gabor23[i][2];
+ descriptor[10]+=this->image->getPixel(vigra::round(x+factor*gabor33[i][0]),vigra::round(y+factor*gabor33[i][1]))*gabor33[i][2];
+ descriptor[11]+=this->image->getPixel(vigra::round(x+factor*gabor43[i][0]),vigra::round(y+factor*gabor43[i][1]))*gabor43[i][2];
+
+ descriptor[12]+=this->image->getPixel(vigra::round(x+factor*gabor53[i][0]),vigra::round(x+factor*gabor53[i][1]))*gabor53[i][2];
+ descriptor[13]+=this->image->getPixel(vigra::round(x+factor*gabor63[i][0]),vigra::round(x+factor*gabor63[i][1]))*gabor63[i][2];
+ descriptor[14]+=this->image->getPixel(vigra::round(x+factor*gabor73[i][0]),vigra::round(x+factor*gabor73[i][1]))*gabor73[i][2];
+ descriptor[15]+=this->image->getPixel(vigra::round(x+factor*gabor83[i][0]),vigra::round(x+factor*gabor83[i][1]))*gabor83[i][2];
+ }
+ for(int i=0; i<17;i++) {
+
+ descriptor[16]+=this->image->getPixel(vigra::round(x+factor*gabor91[i][0]),vigra::round(y+factor*gabor91[i][1]))*gabor91[i][2];
+ descriptor[17]+=this->image->getPixel(vigra::round(x+factor*gabor101[i][0]),vigra::round(y+factor*gabor101[i][1]))*gabor101[i][2];
+ descriptor[18]+=this->image->getPixel(vigra::round(x+factor*gabor111[i][0]),vigra::round(y+factor*gabor111[i][1]))*gabor111[i][2];
+ descriptor[19]+=this->image->getPixel(vigra::round(x+factor*gabor121[i][0]),vigra::round(y+factor*gabor121[i][1]))*gabor121[i][2];
+
+ descriptor[20]+=this->image->getPixel(vigra::round(x+factor*gabor13[i][0]),vigra::round(y+factor*gabor131[i][1]))*gabor131[i][2];
+ descriptor[21]+=this->image->getPixel(vigra::round(x+factor*gabor141[i][0]),vigra::round(y+factor*gabor141[i][1]))*gabor141[i][2];
+ descriptor[22]+=this->image->getPixel(vigra::round(x+factor*gabor151[i][0]),vigra::round(y+factor*gabor151[i][1]))*gabor151[i][2];
+ descriptor[23]+=this->image->getPixel(vigra::round(x+factor*gabor161[i][0]),vigra::round(y+factor*gabor161[i][1]))*gabor161[i][2];
+
+ }
+}
+
+void Descriptor::_RegionResponse(int x,int y, double maxima, double* descriptor) {
+ int add;
+ int sum1;
+ int sum2;
+ int sum3;
+ int sum4;
+ int cnt=0;
+
+ /*for(int i=2; i<=32;i*=2) {
+ add=vigra::round(maxima*i);
+ sum1=this->image->getRegionSum(x-add,y,x,y+add);
+ sum2=-1*this->image->getRegionSum(x,y,x+add,y+add);
+ sum3=this->image->getRegionSum(x,y-add,x+add,y);
+ sum4=-1*this->image->getRegionSum(x-add,y-add,x,y);
+ descriptor[cnt+0]+=(sum1+sum3)/pow((double)i,4.0);
+ descriptor[cnt+1]+=(sum2+sum4)/pow(i,4);
+ descriptor[cnt+2]+=(sum1+sum2+sum3+sum4)/pow(i,4);
+ descriptor[cnt+3]+=abs(sum1+sum2+sum3+sum4)/pow(i,4);
+ cnt+=5;
+ }*/
+}
+
+void Descriptor::_ShapeContext(int x,int y, double maxima, double* descriptor) {
+
+
+ /*for(int i=2; i<=32;i*=2) {
+ add=vigra::round(maxima*i);
+ sum1=this->image->getRegionSum(x-add,y,x,y+add);
+ sum2=-1*this->image->getRegionSum(x,y,x+add,y+add);
+ sum3=this->image->getRegionSum(x,y-add,x+add,y);
+ sum4=-1*this->image->getRegionSum(x-add,y-add,x,y);
+ descriptor[cnt+0]+=(sum1+sum3)/pow((double)i,4.0);
+ descriptor[cnt+1]+=(sum2+sum4)/pow(i,4);
+ descriptor[cnt+2]+=(sum1+sum2+sum3+sum4)/pow(i,4);
+ descriptor[cnt+3]+=abs(sum1+sum2+sum3+sum4)/pow(i,4);
+ cnt+=5;
+ }*/
+}
+
+void Descriptor::createDescriptors() {
+
+// TODO (zoran#8#): Current version detects edges using canny algorithm in the region around the interest point - here is a lot of redundancy. \
+// Need to change: using vigras cannyEdgeImage first run edge detection on the entire picture and then create desciptors.
+ cout << "Creating descriptors..."<< endl;
+ int regionSize;
+ int i=0;
+ double orientation;
+ this->interestPoints->resize(this->interestPoints->size());
+ cout << this->interestPoints->size()<<endl;
+ vector<vector<int> >::iterator iter1 = this->interestPoints->begin();
+ while( iter1 != this->interestPoints->end() ) { //loop over every interest point
+ vector<int > interestPoint=*iter1;
+ i++;
+ vector<double> descriptor;
+ for(int i=0;i<DESCRIPTOR_SIZE;i++) { //initialize array
+ descriptor.push_back(0);
+ }
+
+ regionSize=vigra::round(15*this->_getMaxima(interestPoint[0], interestPoint[1]));
+ double b0,b1;
+ b0=6*this->_getMaxima(interestPoint[0], interestPoint[1]);
+ b1=11*this->_getMaxima(interestPoint[0], interestPoint[1]);
+
+ int xStart=(interestPoint[0]-regionSize);
+ int xEnd=(interestPoint[0] +regionSize);
+ int yStart=(interestPoint[1]-regionSize);
+ int yEnd=(interestPoint[1]+regionSize);
+
+ // only analyse pixels in image
+ int xStartImg = xStart < 0 ? 0 : xStart;
+ int xEndImg = xEnd >= this->image->getWidthBW() ? this->image->getWidthBW() -1 : xEnd;
+ int yStartImg = yStart < 0 ? 0 : yStart;
+ int yEndImg = yEnd >= this->image->getHeightBW() ? this->image->getHeightBW() -1 : yEnd;
+
+/*
+ double result;
+ result = atan (gy/gx) * 180 / PI;
+ result+=90;
+
+ if(result<22.5 && result>=157.5) cp->direction=0;
+ else if(result>=22.5 && result<67.5) cp->direction=1;
+ else if(result>=67.5 && result<112.5) cp->direction=2;
+ else if(result>=112.5 && result<157.5) cp->direction=3;
+
+*/
+ // empty edgel list
+ std::vector<vigra::Edgel> edgels;
+
+ //cout << "Maxima:" << this->_getMaxima(interestPoint[0], interestPoint[1]) << endl;
+ //cout << "Point:" << interestPoint[0] << "," << interestPoint[1]<< endl;
+
+ //
+ vector<int > iPointRelative = interestPoint;
+ iPointRelative[0]=regionSize;
+ iPointRelative[1]=regionSize;
+
+ // find edgels at scale of the interest point
+ vigra::cannyEdgelList(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStartImg, yStartImg),
+ this->image->imgBW->upperLeft() + vigra::Diff2D(xEndImg, yEndImg)),
+ edgels, this->_getMaxima(interestPoint[0], interestPoint[1]),iPointRelative);
+ /*this->image->_cannyEdgelList1(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStart, yStart),
+ this->image->imgBW->upperLeft() + vigra::Diff2D(xEnd, yEnd)),
+ edgels, this->_getMaxima(interestPoint[0], interestPoint[1]), &interestPoint);*/
+ //cout << "Size:" << edgels.size() << endl;
+
+ //first edgel element holds orientation assignment for interest point
+ vector<vigra::Edgel>::iterator iter2 = edgels.begin();
+ vigra::Edgel edgePoint=*iter2;
+ iter2++;
+ //cout << "Orientation:" << edgePoint.orientation << ", strenght:" << edgePoint.strength << endl;
+ orientation = edgePoint.orientation;
+ if(orientation<0) orientation+=2*PI;
+
+ double a0,a1,a2,a3;
+ a0=fmod(0+orientation,2*PI);
+ a1=fmod(PI/2+orientation,2*PI);
+ a2=fmod(PI+orientation,2*PI);
+ a3=fmod((3*PI)/2+orientation,2*PI);
+
+ //int count =0;
+ while( iter2 != edgels.end()) { //loop over every canny pixel
+ edgePoint=*iter2;
+
+ //TODO discard edgels that are not in circle
+ double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x+xStartImg, edgePoint.y+yStartImg);
+
+ if(distance>regionSize) {
+ iter2++; continue;
+ }
+ //count++;
+ //cout << edgePoint.x << ","<<edgePoint.y << "; orientation:"<<edgePoint.orientation<<",strength:"<<edgePoint.strength<<endl;
+ int sector;
+ if(edgePoint.orientation<0) edgePoint.orientation+=2*PI;
+
+ int orientSector;
+
+ //orientation is quantized into 4 bins
+ if(edgePoint.orientation>=0 && edgePoint.orientation<(PI/2)) orientSector=0;
+ else if(edgePoint.orientation>=(PI/2) && edgePoint.orientation<PI) orientSector=1;
+ else if(edgePoint.orientation>=PI && edgePoint.orientation<(3*PI/2)) orientSector=2;
+ else if(edgePoint.orientation>=(3*PI/2)) orientSector=3;
+ else cerr<< "Orientation1:"<<edgePoint.orientation<<endl;
+
+ double result = atan2(regionSize-edgePoint.y,edgePoint.x-regionSize);
+ if(result<0.0) result+=2*PI;
+
+ //location is quantized into 9 bins
+ if(result>=a0 && result<a1) sector=0;
+ else if(result>=a1 && result<a2) sector=1;
+ else if(result>=a2 && result<a3) sector=2;
+ else if(result>=a3 || result < a0) sector=3;
+ else cerr<< "Orientation2:"<<edgePoint.orientation<<endl;
+
+ if(distance<=b0) descriptor[orientSector*9 + 8]+=edgePoint.strength;
+ else if(distance>b0 && distance<=b1) descriptor[orientSector*9 + sector*2+0]+=edgePoint.strength;
+ else if(distance>b1) descriptor[orientSector*9 + sector*2+1]+=edgePoint.strength;
+
+ /*double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x,edgePoint.y);
+ if(distance<=27) descriptor[orientSector*9 + 8]+=edgePoint.strength;
+ else if(distance>27 && distance<=50) descriptor[orientSector*9 + sector*2+0]+=edgePoint.strength;
+ else if(distance>50) descriptor[orientSector*9 + sector*2+1]+=edgePoint.strength;
+*/
+
+ iter2++;
+ }
+ //cout << "Size: "<<edgels.size()<< "; count: "<<count<<endl;
+ this->descriptors.push_back(descriptor);
+ iter1++;
+ }
+}
+vector<vector<double> >* Descriptor::getDescriptors() {
+ return &this->descriptors;
+}
+bool Descriptor::printDescriptors(string name)
+{
+ ofstream o (name.c_str());
+ if (o.is_open()) {
+ o << DESCRIPTOR_SIZE << endl;
+ o << interestPoints->size() << endl;
+ vector<vector<int> >::iterator iter1 = interestPoints->begin();
+ vector<vector<double> >::iterator iterDesc = descriptors.begin();
+
+ //int c=0;
+ while( iter1 != interestPoints->end()) {
+ vector<int > tmp2=*iter1;
+ double r = _getMaxima(tmp2[0], tmp2[1])*HD_INIT_KERNEL_SIZE;
+ o <<tmp2[1]<<" "<<tmp2[0]<<" " << 1/(r*r) << " " << 0 << " " << 1/(r*r) << " ";
+ for (vector<double>::iterator it = iterDesc->begin(); it != iterDesc->end(); ++it) {
+ o << *it << " ";
+ }
+ o << endl;
+ iter1++;
+ iterDesc++;
+ }
+ o.close();
+ return true;
+ } else return false;
+}
+
+bool Descriptor::generateAutopanoXML(string name) {
+ cout << "Generating XML file: "<<name<<endl;
+ ofstream xmlFile (name.c_str());
+ if (xmlFile.is_open())
+ {
+ xmlFile << "<?xml version=\"1.0\" encoding=\"utf-8\"?>"<< endl;
+ xmlFile << "<KeypointXMLList xmlns:xsd=\"http://www.w3.org/2001/XMLSchema\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\">\n";
+ xmlFile << " <XDim>"<< this->image->getWidth() <<"</XDim>"<< endl;
+ xmlFile << " <YDim>"<< this->image->getHeight() <<"</YDim>"<< endl;
+ xmlFile << " <ImageFile>"<< this->image->getPath() <<"</ImageFile>"<< endl;
+ xmlFile << " <Arr>"<< endl;
+
+ vector<vector<int> >::iterator iter1 = interestPoints->begin();
+ vector<vector<double> >::iterator iterDesc = descriptors.begin();
+
+ int c=0;
+ while( iter1 != interestPoints->end()) {
+ vector<int > tmp2=*iter1;
+
+ xmlFile << " <KeypointN>"<< endl;
+ xmlFile << " <X>"<< tmp2[0] <<"</X>"<<endl;
+ xmlFile << " <Y>"<< tmp2[1] <<"</Y>"<<endl;
+ xmlFile << " <Scale>"<< this->_getMaxima(tmp2[0],tmp2[1]) <<"</Scale>"<<endl;
+ xmlFile << " <Orientation>"<< 0 <<"</Orientation>"<<endl;
+ xmlFile << " <Dim>"<< DESCRIPTOR_SIZE <<"</Dim>"<<endl;
+ xmlFile << " <Descriptor>"<<endl;
+
+ for (vector<double>::iterator it = iterDesc->begin(); it != iterDesc->end(); ++it) {
+ xmlFile << " <int>"<< (int)*it << "</int>" << endl;
+ }
+ xmlFile << " </Descriptor>"<<endl;
+ xmlFile << " </KeypointN>"<< endl;
+
+ iter1++;
+ iterDesc++;
+ c++;
+ }
+
+ xmlFile << " </Arr>"<< endl;
+ xmlFile << "</KeypointXMLList>"<< endl;
+
+ xmlFile.close();
+ return true;
+ }
+ else return false;
+}
+
+
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/Descriptor.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,65 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#ifndef DESCRIPTOR_H_INCLUDED
+#define DESCRIPTOR_H_INCLUDED
+#include <string>
+#include "APImage.h"
+#include "HessianDetector.h"
+
+#define NR_ANGLE_BINS 8
+#define DESCRIPTOR_SIZE 36
+
+class Descriptor
+ {
+ public:
+ Descriptor(APImage* i, HessianDetector* hd);
+ double _gaussWeighting(int x, int y, double stdev);
+
+
+
+ void setPoints(vector<vector<int> >* pts);
+ //void orientate();
+ void createDescriptors();
+ vector<vector<double> >* getDescriptors();
+ bool printDescriptors(string name);
+ bool generateAutopanoXML(string name);
+ vector<vector<int> >* interestPoints;
+
+
+ private:
+ APImage* image;
+ HessianDetector* hd;
+
+ double _getMaxima(int x,int y);
+ double _euclidianDistance(int x1, int y1, int x2, int y2);
+ //double _initDescriptor(double* desciptor);
+ void _GaborResponse(int x,int y,double maxima,double* descriptor);
+ void _RegionResponse(int x,int y,double maxima,double* descriptor);
+ void _ShapeContext(int x,int y,double maxima,double* descriptor);
+
+ vector<vector<int> >* maximas;
+ vector<vector<double> > descriptors;
+ //int orientations[];
+ };
+
+
+
+#endif // DESCRIPTOR_H_INCLUDED
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingAlgorithm.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingAlgorithm.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingAlgorithm.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,80 @@
+/**
+ * @file FeatureMatchingAlgorithm.h
+ * gsoc2008_feature_matching
+ *
+ * This class is a superclass for feature matching algorithms
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#ifndef _CTRLPNTSALGORITHMS_FEATURE_MATCHING_ALGO_H
+#define _CTRLPNTSALGORITHMS_FEATURE_MATCHING_ALGO_H
+
+#include <algorithm/PanoramaAlgorithm.h>
+#include <panodata/Panorama.h>
+
+struct FeaturePointMatch
+{
+
+ FeaturePointMatch(unsigned int _im1, unsigned int _im2,
+ Keypoint* _k1, Keypoint* _k2)
+ {
+ im1 = _im1;
+ im2 = _im2;
+ k1 = _k1;
+ k2 = _k2;
+ }
+
+ // variables
+ unsigned int im1, im2;
+ Keypoint k1*, k2*;
+};
+
+// Euclidean distance calculation for two keypoints
+static float eucdist(const Keypoint & p1, const Keypoint & p2)
+{
+ float sum = 0;
+ std::vector<float>::const_iterator it2 = p2.descriptor.begin();
+
+ for(std::vector<float>::const_iterator it1 = p1.descriptor.begin();
+ it1 != p1.descriptor.end(); ++it1, ++it2)
+ {
+ float d = *it1 - *it2;
+ d *= d;
+ sum += d;
+ }
+
+ return pow(sum,0.5f);
+}
+
+// Euclidean distance calculation for two keypoints with an upper bound given
+static float eucdist(const Keypoint & p1, const Keypoint & p2, float ub)
+{
+ float sum = 0;
+ std::vector<float>::const_iterator it2 = p2.descriptor.begin();
+ float ub2 = ub * ub;
+
+ for(std::vector<float>::const_iterator it1 = p1.descriptor.begin();
+ it1 != p1.descriptor.end(); ++it1, ++it2)
+ {
+ float d = *it1 - *it2;
+ d *= d;
+ sum += d;
+
+ // check if the distance exceeds upper bound
+ if (sum > ub2)
+ return numeric_limits<float>::infinity();
+ }
+
+ return pow(sum,0.5f);
+}
+
+class FeatureMatchingAlgorithm : public PanoramaAlgorithm
+{
+
+}
+
+typedef std::vector<FeaturePointMatch> FPMVector;
+
+
+#endif
\ No newline at end of file
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingAlgorithm.h
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,60 @@
+/**
+ * FeatureMatchingLinearSearch.cpp
+ * gsoc2008_feature_matching
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#include "FeatureMatchingLinearSearch.h"
+
+// matching function
+static HuginBase::FPMVector& FeatureMatchingLinearSearch::match(const PanoramaData& pano, float upper_bound_for_distance) {
+
+ FPMVector allMatches;
+ float nearest_keypoint_distance;
+ Keypoint* nearest_point;
+ int nearest_image;
+
+ // for each image in panorama
+ for (unsigned int im=0; im < pano.getNrOfImages(); im++)
+ {
+ // for each keypoint in this image
+ const std::vector<Keypoint> & keypoints = pano.getImage(im).getKeypoints();
+ for (vector<Keypoint>::const_iterator itkey=keypoints.begin();
+ itkey != keypoints.end(); ++itkey)
+ {
+ // find the nearest matching point
+ nearest_keypoint_distance = numeric_limits<float>::infinity();
+ nearest_point = NULL;
+ nearest_image = -1;
+
+ // iterate over the other images
+ for (unsigned int im2=0; im2 < pano.getNrOfImages(); im2++)
+ {
+ if (im2==im) continue;
+ // for each keypoint in this image
+ const std::vector<Keypoint> & keypoints2 = pano.getImage(im2).getKeypoints();
+ for (vector<Keypoint>::const_iterator itkey2=keypoints2.begin();
+ itkey2 != keypoints2.end(); ++itkey2)
+ {
+
+ // calculate euclidean distance
+ float dist = eucdist(*itkey, *itkey2, upper_bound_for_distance);
+ if (dist < nearest_keypoint_distance)
+ {
+ nearest_keypoint_distance = dist;
+ nearest_point = itkey2;
+ nearest_image = im2;
+ }
+ }
+ }
+
+ // if a match found within a neighborhood defined by upperbound,
+ // add to the matching list
+ if (nearest_point != NULL)
+ allMatches.push_back(FeaturePointMatch(im, im2, itkey, itkey2));
+ }
+ }
+
+ return allMatches;
+}
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.cpp
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,53 @@
+/**
+ * @file FeatureMatchingLinearSearch.h
+ * gsoc2008_feature_matching
+ *
+ * This class implements the straightforward way of matching keypoints in a
+ * given panorama.
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#ifndef _CTRLPNTSALGORITHMS_FEATURE_MATCHING_LINEAR_H
+#define _CTRLPNTSALGORITHMS_FEATURE_MATCHING_LINEAR_H
+
+#include <algorithm/PanoramaAlgorithm.h>
+#include <panodata/Panorama.h>
+#include <cmath>
+#include <limits>
+
+class FeatureMatchingLinearSearch : public FeatureMatchingAlgorihtm
+{
+
+public:
+ // constructor and destructor
+ FeatureMatchingLinearSearch(PanoramaData& panorama)
+ : PanoramaAlgorithm(panorama), upper_bound_for_distance(100d) {};
+ FeatureMatchingLinearSearch(PanoramaData& panorama, float ub)
+ : PanoramaAlgorithm(panorama), upper_bound_for_distance(ub) {};
+ virtual ~FeatureMatchingLinearSearch() {};
+
+ // function to find the control points in the panorama
+ virtual bool runAlgorithm()
+ {
+ o_controlPoints = match(o_panorama, upper_bound_for_distance);
+ return true; // let's hope so.
+ }
+
+ // getter method for control points
+ virtual const FPMVector & getControlPoints() const
+ {
+ // [TODO] if(!hasRunSuccessfully()) DEBUG;
+ return o_controlPoints;
+ }
+
+ // matching function
+ static HuginBase::FPMVector& match(const PanoramaData& pano, float upper_bound_for_distance);
+
+protected:
+ FPMVector o_controlPoints;
+ float upper_bound_for_distance;
+};
+
+
+#endif // _H
\ No newline at end of file
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingLinearSearch.h
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,47 @@
+/**
+* FeatureMatchingMultiKdTree.cpp
+ * gsoc2008_feature_matching
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#include "FeatureMatchingMultiKdTree.h"
+
+
+// matching function
+static HuginBase::FPMVector& FeatureMatchingMultiKdTree::match(const PanoramaData& pano) {
+
+ FPMVector allMatches;
+
+ // prepare k-d trees
+ KDTreeKeypointMatcher matcher[pano.getNrOfImages()];
+ for (int im=0; im < pano.getNrOfImages(); im++)
+ {
+ UIntSet single_image;
+ single_image.insert(im);
+ matcher[im] = KDTreeKeypointMatcher(2); // k-d tree with k = 2
+ matcher[im].create(pano, single_image);
+ }
+
+ // keep track of matched points, to avoid re-matching them.
+ UIntSet matchedPoints; // TODO: this strategy cannot be applied to multiple case
+
+ // for each image in panorama
+ for (int im=0; im < pano.getNrOfImages(); im++)
+ {
+ // for each keypoint in this image
+ const std::vector<Keypoint> & keypoints = pano.getImage(im).getKeypoints();
+ for (vector<Keypoint>::const_iterator itkey=keypoints.begin();
+ itkey != keypoints.end(); ++itkey)
+ {
+ // do not try to match this point, if it has already been matched.
+
+ // find all the nearest matching point
+
+ // mark matched points and remove them from further matching.
+
+ }
+ }
+
+ return allMatches;
+}
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.cpp
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,52 @@
+/**
+* @file FeatureMatchingMultiKdTree.h
+ * gsoc2008_feature_matching
+ *
+ * Matching keypoints in a panorama using multiple k-d trees.
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#ifndef _CTRLPNTSALGORITHMS_FEATURE_MATCHING_MULTIKD_H
+#define _CTRLPNTSALGORITHMS_FEATURE_MATCHING_MULTIKD_H
+
+#include <algorithm/PanoramaAlgorithm.h>
+#include <panodata/Panorama.h>
+#include "KDTreeKeypointMatcher.h"
+
+namespace HuginBase {
+
+class FeatureMatchingMultiKdTree : public FeatureMatchingAlgorithm
+{
+
+public:
+ // constructor and destructor
+ FeatureMatchingMultiKdTree(PanoramaData& panorama)
+ : PanoramaAlgorithm(panorama) {};
+ virtual ~FeatureMatchingMultiKdTree() {};
+
+ // function to find the control points in the panorama
+ virtual bool runAlgorithm()
+ {
+ o_controlPoints = match(o_panorama);
+ return true; // let's hope so.
+ }
+
+ // getter method for control points
+ virtual const FPMVector & getControlPoints() const
+ {
+ // [TODO] if(!hasRunSuccessfully()) DEBUG;
+ return o_controlPoints;
+ }
+
+ // matching function
+ static HuginBase::FPMVector& match(const PanoramaData& pano);
+
+protected:
+ FPMVector o_controlPoints;
+};
+
+}
+
+
+#endif
\ No newline at end of file
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingMultiKdTree.h
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,54 @@
+/**
+ * FeatureMatchingSingleKdTree.cpp
+ * gsoc2008_feature_matching
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#include "FeatureMatchingSingleKdTree.h"
+
+namespace HuginBase {
+
+ using namespace std;
+ using namespace hugin_utils;
+
+ // matching function
+ static HuginBase::FPMVector& FeatureMatchingSingleKdTree::match(const PanoramaData& pano) {
+
+ FPMVector allMatches;
+
+ // prepare kd tree
+ KDTreeKeypointMatcher matcher;
+ matcher.create(pano);
+
+ // keep track of matched points, to avoid re-matching them.
+ UIntSet matchedPoints;
+
+ // for each image in panorama
+ for (int im=0; im < pano.getNrOfImages(); im++)
+ {
+ // for each keypoint in this image
+ const std::vector<Keypoint> & keypoints = pano.getImage(*it).getKeypoints();
+ for (vector<Keypoint>::const_iterator itkey=keypoints.begin();
+ itkey != keypoints.end(); ++itkey)
+ {
+ // do not try to match this point, if it has already been matched.
+ if (set_contains(matchedPoints, pnr)) continue;
+
+ // find the nearest matching point
+ ImageKeypoint onematch = matcher.match(*itkey, *it);
+ if (onematch != NULL)
+ allMatches.push_back(FeaturePointMatch(im, onematch.imageNr, itkey, onematch.keypoint));
+
+ // mark matched points and remove them from further matching.
+ for (unsigned j=0; j < matches.size(); j++) {
+ // mark controlpoint as matched.
+ matchedPoints.insert(matcher.getKeypointIdxOfMatch(j));
+ }
+ }
+ }
+
+ return allMatches;
+ }
+
+}
\ No newline at end of file
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.cpp
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,52 @@
+/**
+ * @file FeatureMatchingSingleKdTree.h
+ * gsoc2008_feature_matching
+ *
+ * Matching keypoints in a panorama using single k-d tree.
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#ifndef _CTRLPNTSALGORITHMS_FEATURE_MATCHING_SINGLEKD_H
+#define _CTRLPNTSALGORITHMS_FEATURE_MATCHING_SINGLEKD_H
+
+#include <algorithm/PanoramaAlgorithm.h>
+#include <panodata/Panorama.h>
+#include "KDTreeKeypointMatcher.h"
+
+namespace HuginBase {
+
+ class FeatureMatchingSingleKdTree : public FeatureMatchingAlgorithm
+ {
+
+ public:
+ // constructor and destructor
+ FeatureMatchingSingleKdTree(PanoramaData& panorama)
+ : PanoramaAlgorithm(panorama) {};
+ virtual ~FeatureMatchingSingleKdTree() {};
+
+ // function to find the control points in the panorama
+ virtual bool runAlgorithm()
+ {
+ o_controlPoints = match(o_panorama);
+ return true; // let's hope so.
+ }
+
+ // getter method for control points
+ virtual const FPMVector & getControlPoints() const
+ {
+ // [TODO] if(!hasRunSuccessfully()) DEBUG;
+ return o_controlPoints;
+ }
+
+ // matching function
+ static HuginBase::FPMVector& match(const PanoramaData& pano);
+
+ protected:
+ FPMVector o_controlPoints;
+ };
+
+}
+
+
+#endif
\ No newline at end of file
Property changes on: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/FeatureMatchingSingleKdTree.h
___________________________________________________________________
Name: svn:executable
+ *
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,695 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <limits>
+#include "HessianDetector.h"
+
+using namespace std;
+
+HessianDetector::HessianDetector(APImage* i, int nrPoints, CONVOLUTION_TYPE type, int octaves) {
+ if(octaves>HD_MAX_OCTAVES) return; // TODO (zoran#1#): Throw error or log status.
+
+ this->image=i;
+ this->nrPoints=nrPoints;
+ this->nrOctaves=octaves;
+ this->convolutionType = type;
+
+ vector<int> point;
+ point.resize(2);
+
+ this->orderedList.reserve(sizeof(point)*nrPoints);
+
+ //initialize the vectors
+ determinants.resize(this->image->getWidthBW());
+ maximas.resize(this->image->getWidthBW());
+
+ const int iMax=std::numeric_limits<int>::max();
+ for(int i=0;i<this->image->getWidthBW();i++) {
+ maximas[i].resize(this->image->getHeightBW());
+ determinants[i].resize(this->image->getHeightBW());
+ for(int j=0; j<this->image->getHeightBW();j++) {
+ maximas[i][j] =0 ;
+ determinants[i][j]=iMax;
+ }
+ }
+}
+
+bool HessianDetector::detect() {
+ if(this->convolutionType==HD_BOX_FILTERS) {
+ return this->_boxFilterDetect();
+ } else if(this->convolutionType==HD_SLIDING_WINDOW) {
+ return this->_slidingWDetect();
+ }
+ return false;
+}
+
+bool HessianDetector::_slidingWDetect() {
+ int height = this->image->getHeightBW();
+ int width = this->image->getWidthBW();
+
+ int kernelX=9;
+ int kernelY=9;
+ int kernelSize=9;
+ double scale=1.2;
+
+ int gxx[][9]={
+ 1, 7, 29, 68, 90, 68, 29, 7, 1, //9x9 gauss kernels
+ 4, 26, 106, 247, 327, 247, 106, 26, 4,
+ 5, 33, 133, 310, 411, 310, 133, 33, 5,
+ -4, -27, -109, -252, -335, -252, -109, -27, -4,
+ -11, -81, -329, -765, -1013, -765, -329, -81, -11,
+ -4, -27, -109, -252, -335, -252, -109, -27, -4,
+ 5, 33, 133, 310, 411, 310, 133, 33, 5,
+ 4, 26, 106, 247, 327, 247, 106, 26, 4,
+ 1, 7, 29, 68, 90, 68, 29, 7, 1
+ };
+
+ int gyy[][9]={
+ 1, 4, 5, -4, -11, -4, 5, 4, 1,
+ 7, 26, 33, -27, -81, -27, 33, 26, 7,
+ 29, 106, 133, -109, -329, -109, 133, 106, 29,
+ 68, 247, 310, -252, -765, -252, 310, 247, 68,
+ 90, 327, 411, -335, -1013, -335, 411, 327, 90,
+ 68, 247, 310, -252, -765, -252, 310, 247, 68,
+ 29, 106, 133, -109, -329, -109, 133, 106, 29,
+ 7, 26, 33, -27, -81, -27, 33, 26, 7,
+ 1, 4, 5, -4, -11, -4, 5, 4, 1
+ };
+
+ int gxy[][9]={
+ 1, 5, 15, 17, 0, -17, -15, -5, -1,
+ 5, 29, 78, 91, 0, -91, -78, -29, -5,
+ 15, 78, 214, 248, 0, -248, -214, -78, -15,
+ 17, 91, 248, 289, 0, -289, -248, -91, -17,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ -17, -91, -248, -289, 0, 289, 248, 91, 17,
+ -15, -78, -214, -248, 0, 248, 214, 78, 15,
+ -5, -29, -78, -91, 0, 91, 78, 29, 5,
+ -1, -5, -15, -17, 0, 17, 15, 5, 1
+ };
+
+ int offsetX=floor((double)kernelX/2);
+ int offsetY=floor((double)kernelY/2);
+ int xStart=0;
+ int yStart=0;
+ int pixelSumYY;
+ int pixelSumXX;
+ int pixelSumXY;
+
+ int tmpX=0;
+ int tmpY=0;
+ int val;
+ int mappingX;
+ int mappingY;
+
+ int det;
+
+
+ for(int s=0; s<6;s+=1.0) { //scale-space loop
+ //cout << s << "\n";
+ height=image->getHeightBW();
+ width=image->getWidthBW();
+ /*cout << height << "," << width << " \n";
+ cout << " ";*/
+ for(int i=0;i<height;i++) {
+ for(int j=0; j<width;j++) {
+ pixelSumYY=0;
+ pixelSumXX=0;
+ pixelSumXY=0;
+
+ xStart=i-offsetX;
+ yStart=j-offsetY;
+
+ for(int k=0;k<kernelX;k++) {
+ tmpX=xStart+k;
+ if(tmpX<0) tmpX=height+tmpX;
+ if(tmpX>=height) tmpX=tmpX-height;
+
+ for(int l=0;l<kernelY;l++) {
+ tmpY=yStart+l;
+
+ if(tmpY<0) tmpY=width+tmpY;
+ if(tmpY>=height) tmpY=tmpY-height;
+
+ val=image->getPixel(tmpX,tmpY);
+
+ pixelSumYY+=gyy[k][l]*val;
+ pixelSumXX+=gxx[k][l]*val;
+ pixelSumXY+=gxy[k][l]*val;
+ }
+ }
+
+ //determinant of the Hessian Matrix
+ det = pixelSumXX*pixelSumYY-pow((double)pixelSumXY,2.0);
+
+ //map coordinates on the scaled image onto coordinates on the original image
+ mappingX=vigra::round(i*pow(scale,s));
+ mappingY=vigra::round(j*pow(scale,s));
+
+ if(det>determinants[mappingX][mappingY]) {
+ //if(det>max) max=det;
+ determinants[mappingX][mappingY]=det;
+ maximas[mappingX][mappingY]=kernelSize*pow(scale,s);
+ }
+ }
+ //cout << "\n";
+ }
+ image->scale(scale);
+ } //scale-space
+
+ int xEnd=0;
+ int yEnd=0;
+
+ for(int i=0;i<image->getHeight();i++) {
+ if(i>=(image->getHeight()-1)) {
+ xStart=i-1;
+ xEnd=0;
+ } else if(i==0) {
+ xStart=image->getHeight()-1;
+ xEnd=i+1;
+ } else {
+ xStart=i-1;
+ xEnd=i+1;
+ }
+ for(int j=0; j<image->getWidth();j++) {
+ if(j>=(image->getWidth()-1)) {
+ yStart=j-1;
+ yEnd=0;
+ } else if(j==0) {
+ yStart=image->getWidth()-1;
+ yEnd=j+1;
+ } else {
+ yStart=j-1;
+ yEnd=j+1;
+ }
+ /**if any determinant in the area around the pixel is
+ *greater than the pixel, suppress the current pixel
+ *
+ * This can also be done in a loop.
+ *
+ */
+
+ if(determinants[i][j]>determinants[xStart][yStart] &&
+ determinants[i][j]>determinants[xStart][j] &&
+ determinants[i][j]>determinants[xStart][yEnd] &&
+ determinants[i][j]>determinants[i][yStart] &&
+ determinants[i][j]>determinants[i][yEnd] &&
+ determinants[i][j]>determinants[xEnd][yStart] &&
+ determinants[i][j]>determinants[xEnd][j] &&
+ determinants[i][j]>determinants[xEnd][yEnd] ) {
+ this->_insertToList(&i,&j);
+ }
+
+ }
+ //cout << "\n";
+ }
+ return true;
+}
+
+bool HessianDetector::_boxFilterDetect() {
+
+ int height = this->image->getHeightBW();
+ int width = this->image->getWidthBW();
+
+ cout << "Height:"<<height << ",width:"<<width<<endl;
+ //cout << "Max vector size:"<< orderedList.max_size() << ", capacity:" << orderedList.capacity() << "\n";
+
+
+ const int iMax=std::numeric_limits<int>::max();
+
+ int offsetX=width*0.05;
+ int offsetY=height*0.05;
+
+
+ int xStart=offsetX;
+ int yStart=offsetY;
+ int xEnd=width-offsetX;
+ int yEnd=height-offsetY;
+
+ for(int j=yStart;j<yEnd;j++) {
+ for(int i=xStart; i<xEnd;i++) {
+ //_calculateMaxDet(x_coord,y_coord)
+ _calculateMaxDet(i,j); //calculate determinant
+
+ }
+ }
+
+ vector<vector<int> >::iterator it;
+
+ it = orderedList.begin();
+ double sum=0;
+ double count=0;
+ //point (x,y) => (j,i)
+ for(int j=yStart;j<yEnd;j++) {
+ for(int i=xStart; i<xEnd;i++) {
+ /**Non-maxima suppression
+ *if any determinant in the area around the pixel is
+ *greater than the pixel, suppress the current pixel.
+ *
+ * This can also be done in a loop.
+ *
+ */
+
+ if(determinants[i][j]>determinants[i-1][j-1] &&
+ determinants[i][j]>determinants[i-1][j] &&
+ determinants[i][j]>determinants[i-1][j+1] &&
+ determinants[i][j]>determinants[i][j-1] &&
+ determinants[i][j]>determinants[i][j+1] &&
+ determinants[i][j]>determinants[i+1][j-1] &&
+ determinants[i][j]>determinants[i+1][j] &&
+ determinants[i][j]>determinants[i+1][j+1] &&
+ determinants[i][j]>300000) {
+ //this->_insertToList(&i,&j);
+ //cout << maximas[i][j]<<","<<determinants[i][j]<< ";"<<endl;
+ /*cout << i<<","<<j<< "\n";*/
+ vector<int> point;
+ point.resize(2);
+
+ point[0]=i;
+ point[1]=j;
+
+ this->orderedList.push_back(point);
+
+ count++;
+ sum+=determinants[i][j];
+ //it = orderedList.insert ( it , point);
+
+ }
+ //cout << maximas[i][j]<<","<<determinants[i][j]<< " ";
+ }
+ //cout << "\n";
+ }
+ cout << "Detected points:"<< this->orderedList.size() <<" ";
+ cout << "Average:"<< sum/count <<endl;
+ this->_cutPointList(sum/count,10000);
+ return true;
+}
+
+void HessianDetector::_cutPointList(double threshold, int nrPoints) {
+ if(this->orderedList.size()<nrPoints) return;
+ vector<vector<int> >::iterator iter1 = this->orderedList.begin();
+ double avg=0;
+ int count=0;
+ while( iter1 != this->orderedList.end()) { //loop over every interest point
+ vector<int > interestPoint=*iter1;
+
+ iter1++;
+ vector<int > interestPoint2=*iter1;
+
+ //cout << this->determinants[interestPoint[0]][interestPoint[1]] << endl;
+ //cout << count << endl;
+
+
+ iter1--;
+
+ if(this->determinants[interestPoint[0]][interestPoint[1]]<threshold) {
+ this->orderedList.erase(iter1);
+ } else {
+ avg+=this->determinants[interestPoint[0]][interestPoint[1]];
+ iter1++;
+ }
+ count++;
+ }
+
+ cout << "Number of points:"<< this->orderedList.size()<< endl;
+
+ this->_cutPointList(avg/this->orderedList.size(), nrPoints);
+}
+
+void HessianDetector::_calculateMaxDet(int x, int y) {
+ /**
+ * Each octave has 4 kernel sizes.
+ */
+ int octaves[][4] ={
+ 9,15,21,27,
+ 15,27,39,51,
+ 21,45,69,93
+ };
+ double interp1, interp2;
+
+ //first octave
+ interp1=(this->_convolutePixel(&x,&y,&octaves[0][0])+this->_convolutePixel(&x,&y,&octaves[0][1]))/2;
+ interp2=(this->_convolutePixel(&x,&y,&octaves[0][2])+this->_convolutePixel(&x,&y,&octaves[0][3]))/2;
+ if(interp1>interp2) {
+ this->determinants[x][y]=interp1;
+ this->maximas[x][y]=(((octaves[0][0]+octaves[0][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
+ } else {
+ this->determinants[x][y]=interp2;
+ this->maximas[x][y]=(((octaves[0][2]+octaves[0][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
+ }
+
+ //second octave
+ if(this->nrOctaves>1 && x%2==0 && y%2==0) { //points are sampled
+ interp1=(this->_convolutePixel(&x,&y,&octaves[1][0])+this->_convolutePixel(&x,&y,&octaves[1][1]))/2;
+ interp2=(this->_convolutePixel(&x,&y,&octaves[1][2])+this->_convolutePixel(&x,&y,&octaves[1][3]))/2;
+ if(interp1>this->determinants[x][y] || interp2>this->determinants[x][y] ) {
+ if(interp1>interp2) {
+ this->determinants[x][y]=interp1;
+ this->maximas[x][y]=(((octaves[1][0]+octaves[1][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
+ } else {
+ this->determinants[x][y]=interp2;
+ this->maximas[x][y]=(((octaves[1][2]+octaves[1][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
+ }
+ }
+ }
+
+ //third octave
+ if(this->nrOctaves==HD_MAX_OCTAVES && x%4==0 && y%4==0) { //points are sampled
+ interp1=(this->_convolutePixel(&x,&y,&octaves[2][0])+this->_convolutePixel(&x,&y,&octaves[2][1]))/2;
+ interp2=(this->_convolutePixel(&x,&y,&octaves[2][2])+this->_convolutePixel(&x,&y,&octaves[2][3]))/2;
+ if(interp1>this->determinants[x][y] || interp2>this->determinants[x][y] ) {
+ if(interp1>interp2) {
+ this->determinants[x][y]=interp1;
+ this->maximas[x][y]=(((octaves[2][0]+octaves[2][1])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE;
+ } else {
+ this->determinants[x][y]=interp2;
+ this->maximas[x][y]=(((octaves[2][2]+octaves[2][3])/2)*HD_INITIAL_SCALE)/HD_INIT_KERNEL_SIZE; //scale of interpolation
+ }
+ }
+ }
+}
+
+int HessianDetector::_getHessianDeterminant(int* pixelSumXX, int* pixelSumXY, int* pixelSumYY) {
+ return *pixelSumXX*(*pixelSumYY)-pow((0.9*(*pixelSumXY)),2);
+}
+
+int HessianDetector::_convolutePixel(int* x, int* y, int* kernelSize) {
+//for box filters only
+ int pixelSumXX=0;
+ int pixelSumYY=0;
+ int pixelSumXY=0;
+ int det;
+ switch(*kernelSize) {
+ case 9:
+ //9x9 kernel
+ //getRegionSum(y1,x1,y2,x2)
+ pixelSumYY=this->image->getRegionSum(*y-4,*x-2,*y-2,*x+2);
+ pixelSumYY+=this->image->getRegionSum(*y+2,*x-2,*y+4,*x+2);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-1,*x-2,*y+1,*x+2);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-3,*x+1,*y-1,*x+3);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-3,*y+3,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-3,*x-3,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+3,*x+3);
+
+ pixelSumXX=this->image->getRegionSum(*y-2,*x-4,*y+2,*x-2);
+ pixelSumXX+=this->image->getRegionSum(*y-2,*x+2,*y+2,*x+4);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-2,*x-1,*y+2,*x+1);
+
+ //determinant of the Hessian Matrix
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY);
+ break;
+ case 11:
+ //11x11 kernel
+ pixelSumYY=this->image->getRegionSum(*y-5,*x-3,*y-2,*x+3);
+ pixelSumYY+=this->image->getRegionSum(*y+2,*x-3,*y+5,*x+3);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-1,*x-3,*y+1,*x+3);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-4,*x+1,*y-1,*x+4);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-4,*y+4,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-4,*x-4,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+4,*x+4);
+
+ pixelSumXX=this->image->getRegionSum(*y-3,*x-5,*y+3,*x-2);
+ pixelSumXX+=this->image->getRegionSum(*y-3,*x+2,*y+3,*x+5);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-3,*x-1,*y+3,*x+1);
+ //determinant of the Hessian Matrix
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/2.232;
+ break;
+ case 15:
+ //15x15 kernel
+ pixelSumYY=this->image->getRegionSum(*y-7,*x-4,*y-3,*x+4);
+ pixelSumYY+=this->image->getRegionSum(*y+3,*x-4,*y+7,*x+4);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-2,*x-4,*y+2,*x+4);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-5,*x+1,*y-1,*x+5);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-5,*y+5,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-5,*x-5,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+5,*x+5);
+
+ pixelSumXX=this->image->getRegionSum(*y-4,*x-7,*y+4,*x-3);
+ pixelSumXX+=this->image->getRegionSum(*y-4,*x+3,*y+4,*x+7);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-4,*x-2,*y+4,*x+2);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/7.72;
+ break;
+ case 17:
+ //17x17 kernel
+ pixelSumYY=this->image->getRegionSum(*y-8,*x-4,*y-3,*x+4);
+ pixelSumYY+=this->image->getRegionSum(*y+3,*x-4,*y+8,*x+4);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-2,*x-4,*y+2,*x+4);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-6,*x+1,*y-1,*x+6);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-6,*y+6,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-6,*x-6,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+6,*x+6);
+
+ pixelSumXX=this->image->getRegionSum(*y-4,*x-8,*y+4,*x-3);
+ pixelSumXX+=this->image->getRegionSum(*y-4,*x+3,*y+4,*x+8);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-4,*x-2,*y+4,*x+2);
+ //determinant of the Hessian Matrix
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/12.73;
+ break;
+ case 21:
+ //21x21 kernel
+ pixelSumYY=this->image->getRegionSum(*y-10,*x-6,*y-4,*x+6);
+ pixelSumYY+=this->image->getRegionSum(*y+4,*x-6,*y+10,*x+6);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-3,*x-6,*y+3,*x+6);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-7,*x+1,*y-1,*x+7);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-7,*y+7,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-7,*x-7,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+7,*x+7);
+
+ pixelSumXX=this->image->getRegionSum(*y-6,*x-10,*y+6,*x-4);
+ pixelSumXX+=this->image->getRegionSum(*y-6,*x+4,*y+6,*x+10);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-6,*x-3,*y+6,*x+3);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/29.642;
+ break;
+ case 27:
+ //27x27 kernel
+ pixelSumYY=this->image->getRegionSum(*y-13,*x-7,*y-5,*x+7);
+ pixelSumYY+=this->image->getRegionSum(*y+5,*x-7,*y+13,*x+7);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-4,*x-7,*y+4,*x+7);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-9,*x+1,*y-1,*x+9);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-9,*y+9,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-9,*x-9,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+9,*x+9);
+
+ pixelSumXX=this->image->getRegionSum(*y-7,*x-13,*y+7,*x-5);
+ pixelSumXX+=this->image->getRegionSum(*y-7,*x+5,*y+7,*x+13);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-7,*x-4,*y+7,*x+4);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/81;
+ break;
+ case 39:
+ //39x39 kernel
+ pixelSumYY=this->image->getRegionSum(*y-19,*x-11,*y-7,*x+11);
+ pixelSumYY+=this->image->getRegionSum(*y+7,*x-11,*y+19,*x+11);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-6,*x-11,*y+6,*x+11);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-13,*x+1,*y-1,*x+13);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-13,*y+13,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-13,*x-13,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+13,*x+13);
+
+ pixelSumXX=this->image->getRegionSum(*y-11,*x-19,*y+11,*x-7);
+ pixelSumXX+=this->image->getRegionSum(*y-11,*x+7,*y+11,*x+19);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-11,*x-6,*y+11,*x+6);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/352.61;
+ break;
+ case 45:
+ //45x45 kernel
+ pixelSumYY=this->image->getRegionSum(*y-22,*x-12,*y-8,*x+12);
+ pixelSumYY+=this->image->getRegionSum(*y+8,*x-12,*y+22,*x+12);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-7,*x-12,*y+7,*x+12);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-15,*x+1,*y-1,*x+15);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-15,*y+15,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-15,*x-15,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+15,*x+15);
+
+ pixelSumXX=this->image->getRegionSum(*y-12,*x-22,*y+12,*x-8);
+ pixelSumXX+=this->image->getRegionSum(*y-12,*x+8,*y+12,*x+22);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-12,*x-7,*y+12,*x+7);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/625;
+ break;
+ case 51:
+ //51x51 kernel
+ pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
+ pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
+
+ pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
+ pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/1031.123;
+ break;
+ case 69:
+// TODO (zoran#1#):
+ //69x69 kernel
+ pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
+ pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
+
+ pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
+ pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/32.1;
+ break;
+ case 93:
+// TODO (zoran#1#):
+ //93x93 kernel
+ pixelSumYY=this->image->getRegionSum(*y-25,*x-14,*y-9,*x+14);
+ pixelSumYY+=this->image->getRegionSum(*y+9,*x-14,*y+25,*x+14);
+ pixelSumYY+=-2*this->image->getRegionSum(*y-8,*x-14,*y+8,*x+14);
+
+ pixelSumXY=-1*this->image->getRegionSum(*y-17,*x+1,*y-1,*x+17);
+ pixelSumXY+=-1*this->image->getRegionSum(*y+1,*x-17,*y+17,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y-17,*x-17,*y-1,*x-1);
+ pixelSumXY+=this->image->getRegionSum(*y+1,*x+1,*y+17,*x+17);
+
+ pixelSumXX=this->image->getRegionSum(*y-14,*x-25,*y+14,*x-9);
+ pixelSumXX+=this->image->getRegionSum(*y-14,*x+9,*y+14,*x+25);
+ pixelSumXX+=-2*this->image->getRegionSum(*y-14,*x-8,*y+14,*x+8);
+
+ det=_getHessianDeterminant(&pixelSumXX,&pixelSumXY,&pixelSumYY)/32.1;
+ break;
+ }
+ return det;
+}
+
+void HessianDetector::_insertToList(int* x, int* y) {
+ vector<int> point;
+ point.resize(2);
+ point[0]=*x;
+ point[1]=*y;
+
+ if(this->orderedList.size()>=this->nrPoints) {
+ vector<int> tmp=this->orderedList.back();
+ if(determinants[*x][*y]<=determinants[tmp[0]][tmp[1]]) return;
+
+ tmp = this->orderedList.front();
+ if(determinants[*x][*y]>=determinants[tmp[0]][tmp[1]]) {
+ this->orderedList.insert(this->orderedList.begin(),point);
+ this->orderedList.pop_back();
+ return;
+ }
+ this->orderedList.pop_back();
+ } else {
+ if(this->orderedList.size()==0) {
+ this->orderedList.push_back(point);
+ return;
+ }
+ vector<int> tmp=this->orderedList.back();
+ if(determinants[*x][*y]<=determinants[tmp[0]][tmp[1]]) {
+ this->orderedList.push_back(point);
+ return;
+ }
+ tmp = this->orderedList.front();
+ if(determinants[*x][*y]>=determinants[tmp[0]][tmp[1]]) {
+ orderedList.insert(orderedList.begin(),point);
+ return;
+ }
+ }
+ vector<vector<int> >::iterator iter1 = orderedList.begin();
+ while( iter1 != orderedList.end()) {
+ vector<int > tmp=*iter1;
+ if(determinants[tmp[0]][tmp[1]]<=determinants[*x][*y]) {
+ orderedList.insert(iter1,point);
+ return;
+ }
+ iter1++;
+ }
+}
+void HessianDetector::printPoints() {
+ vector<vector<int> >::iterator iter1 = orderedList.begin();
+ //int c=0;
+ while( iter1 != orderedList.end()) {
+ vector<int > tmp2=*iter1;
+ cout << "("<<tmp2[0]<<","<<tmp2[1]<<","<<tmp2[2]<<")\n";
+ this->image->drawCircle(tmp2[0],tmp2[1],maximas[tmp2[0]][tmp2[1]]*HD_INIT_KERNEL_SIZE);
+ //c++;
+ //this->image->drawCircle(tmp2[1],tmp2[0],1);
+ iter1++;
+ }
+ this->image->show();
+ //cout << c <<"\n";
+}
+void HessianDetector::printPoints(std::ostream & o) {
+ cout << "Print HD points"<<endl;
+
+ vector<vector<int> >::iterator iter1 = orderedList.begin();
+ int n=0;
+
+/*
+ while( iter1 != orderedList.end()) {
+ n++;
+ iter1++;
+ }*/
+ o << 1 << endl;
+ o << this->orderedList.size() << endl;
+ iter1 = orderedList.begin();
+ //int c=0;
+ while( iter1 != orderedList.end()) {
+ vector<int > tmp2=*iter1;
+ //double r = getMaxima(tmp2[0], tmp2[1])*HD_INIT_KERNEL_SIZE;
+ double r = 1;
+ o <<tmp2[1]<<" "<<tmp2[0]<<" "<< " " << 1/(r*r) << " " << 0 << " " << 1/(r*r) << endl;
+ iter1++;
+ }
+ cout << "EOF Print HD points"<<endl;
+
+}
+
+vector<vector<int> >* HessianDetector::getPoints() {
+ return &this->orderedList;
+}
+int HessianDetector::getNrPoints() {
+ return this->nrPoints;
+}
+double HessianDetector::getMaxima(int x, int y) {
+ return this->maximas[x][y];
+}
+double HessianDetector::_getScale(int kernelSize) {
+ return kernelSize/HD_INIT_KERNEL_SIZE;
+}
+void HessianDetector::dump() {
+ delete(&determinants);
+}
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.h
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.h (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/HessianDetector.h 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,86 @@
+/***************************************************************************
+ * Copyright (C) 2007 by Zoran Mesec *
+ * zoran.mesec@... *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#ifndef HESSIANDETECTOR_H_
+#define HESSIANDETECTOR_H_
+#include "APImage.h"
+
+using namespace std;
+
+#define CONVOLUTION_TYPE int
+#define HD_SLIDING_WINDOW 2
+#define HD_BOX_FILTERS 1
+#define HD_MAX_OCTAVES 3
+#define HD_INITIAL_SCALE 1.2
+
+//initial kernel size, representing scale 1.2
+#define HD_INIT_KERNEL_SIZE 9
+
+#define PI 3.14159265
+
+class HessianDetector
+ {
+ public:
+ HessianDetector(APImage* i, int nrPoints=1000, CONVOLUTION_TYPE type=HD_BOX_FILTERS, int nrOctaves=1);
+
+ bool detect();
+ void printPoints();
+ void printPoints(std::ostream & o);
+
+ //clears memory of data structures that are not required anymore(before the description process)
+ void dump();
+ vector<vector<int> >* getPoints();
+ double getMaxima(int x, int y);
+ int getNrPoints();
+
+/* private slots:
+*/
+
+ private:
+ APImage* image;
+ int nrPoints;
+ int nrOctaves;
+ CONVOLUTION_TYPE convolutionType;
+
+ /**
+ * 0---> x
+ * |
+ * |
+ * v y
+ */
+ //Point is represented as T(x,y)
+ //determinants[x][y] stands for determinant of point T(x,y)
+ std::vector<vector<int> > determinants; //holds the values of the hessian determinant for each pixel
+ std::vector<vector<int> > orderedList; //first n pixels with the largest determinant values
+ //determinants[x][y]
+ std::vector<vector<double> > maximas; //holds the scales where scale-space representation for each pixel attends maximum
+
+ int _getHessianDeterminant(int* pixelSumXX, int* pixelSumXY, int *pixelSumYY);
+ void _calculateMaxDet(int i, int j); //calculates scale-space maxima for pixel at coord i,j
+ void _cutPointList(double average, int nrPoints);
+ double _getScale(int kernelSize);
+ int _convolutePixel(int* coordX, int* coordY, int* kernelSize);
+
+ bool _boxFilterDetect();
+ bool _slidingWDetect();
+ void _insertToList(int* x, int* y);
+ };
+
+#endif /*HESSIANDETECTOR_H_*/
Added: hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/KDTreeKeypointMatcher.cpp
===================================================================
--- hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/KDTreeKeypointMatcher.cpp (rev 0)
+++ hugin/branches/gsoc2008_feature_matching/src/hugin_base/algorithms/control_points/KDTreeKeypointMatcher.cpp 2008-06-28 11:33:30 UTC (rev 3154)
@@ -0,0 +1,143 @@
+/**
+ * @file KDTreeKeypointMatcher.cpp
+ * gsoc2008_feature_matching
+ *
+ * This class is a generic k-d tree wrapper
+ *
+ * @author Onur Kucuktunc <onurcc@...>
+ */
+
+#include "KDTreeKeypointMatcher.h"
+
+/**
+* create KDTree from the feature descriptors stored in the panorama.
+*
+* @param pano Panorama from which the keypoints are taken.
+*/
+void KDTreeKeypointMatcher::create(const PanoramaData & pano)
+{
+ UIntSet imgs;
+ // for each image in panorama
+ for (int im=0; im < pano.getNrOfImages(); im++)
+ imgs.insert(im);
+
+ // call the create function with all images
+ create(pano, &imgs);
+}
+
+/**
+* create KDTree from the feature descriptors stored in the panorama.
+*
+* @param pano Panorama from which the keypoints are taken.
+* @param images Only use keypoints from the specified images.
+*/
+void KDTreeKeypointMatcher::create(const PanoramaData & pano, const UIntSet & imgs)
+{
+ int nKeypoints=0;
+ int dim=0;
+
+ // first, count how many descriptors we have, and the dimension of the space
+ for (UIntSet::const_iterator it=imgs.begin(); it != imgs.end(); ++it) {
+ const std::vector<Keypoint> & keypoints = pano.getImage(*it).getKeypoints();
+ nKeypoints += keypoints.size();
+ if (keypoints.size() > 0)
+ dim = keypoints[0].descriptor.size();
+ }
+ assert(dim > 0);
+
+ // prepare a big list of points, and insert all points
+ if (m_allPoints) delete[] m_allPoints;
+ m_allPoints = new ANNpoint[nKeypoints];
+
+ ANNpointArray pointsPtr = m_allPoints;
+
+ // for all images
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|