summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMichael M <myrhev@Myrhev-laptop.(none)>2009-11-19 13:40:29 +0100
committerMichael M <myrhev@Myrhev-laptop.(none)>2009-11-19 13:40:29 +0100
commitc503c0dd53d74c8807be0a75e0c0e4ccec77a935 (patch)
tree6c3eae2c30140575a3dfe345a70c17faebe90bec
downloadpanor-c503c0dd53d74c8807be0a75e0c0e4ccec77a935.tar.gz
panor-c503c0dd53d74c8807be0a75e0c0e4ccec77a935.tar.bz2
Initial add
-rw-r--r--CMakeLists.txt10
-rw-r--r--Imagine/Features.h28
-rw-r--r--Imagine/Features/Feat.h54
-rw-r--r--Imagine/Features/IO.h209
-rw-r--r--Imagine/Features/Match.h106
-rw-r--r--Imagine/Features/SIFT.h77
-rw-r--r--Imagine/Optim.h20
-rw-r--r--Imagine/Optim/Ransac.h207
-rw-r--r--display.cpp37
-rw-r--r--display.h14
-rw-r--r--image0006.jpgbin0 -> 177671 bytes
-rw-r--r--image0007.jpgbin0 -> 184444 bytes
-rw-r--r--image0008.jpgbin0 -> 200337 bytes
-rw-r--r--images.cpp17
-rw-r--r--images.h33
-rw-r--r--main.cpp54
-rw-r--r--matches.cpp28
-rw-r--r--matches.h21
-rw-r--r--vl/generic.c400
-rw-r--r--vl/generic.h255
-rw-r--r--vl/host.c509
-rw-r--r--vl/host.h367
-rw-r--r--vl/imop.c85
-rw-r--r--vl/imop.h62
-rw-r--r--vl/imop.tc155
-rw-r--r--vl/mathop.c30
-rw-r--r--vl/mathop.h450
-rw-r--r--vl/mser.c971
-rw-r--r--vl/mser.h422
-rw-r--r--vl/sift.c2175
-rw-r--r--vl/sift.h410
31 files changed, 7206 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..1dc071f
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,10 @@
+PROJECT(Panor)
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
+SET(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} $ENV{IMAGINEPP_ROOT}/CMake)
+FIND_PACKAGE(Imagine)
+
+SET(FEATS SIFT.cpp vl/imop.c vl/sift.c vl/generic.c vl/host.c)
+
+ImagineAddExecutable(Matches main.cpp display.cpp images.cpp matches.cpp framework.cpp ${FEATS})
+ImagineUseModules(Matches Images LinAlg)
diff --git a/Imagine/Features.h b/Imagine/Features.h
new file mode 100644
index 0000000..0bd9b33
--- /dev/null
+++ b/Imagine/Features.h
@@ -0,0 +1,28 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+#ifndef _IMAGINEFEATURES_H
+#define _IMAGINEFEATURES_H
+
+#include <string>
+#include <fstream>
+#include <list>
+
+#include <Imagine/Images.h>
+
+/// \defgroup Features Features Library
+/// @{
+/// \example Features/test/test.cpp
+/// \example Features/test/example.cpp
+/// @}
+
+#include "Features/Feat.h" // Main class
+#include "Features/SIFT.h" // VLFeat SIFT
+#include "Features/Match.h" // Basic Matching
+#include "Features/IO.h" // Files + Display
+
+#endif
+
diff --git a/Imagine/Features/Feat.h b/Imagine/Features/Feat.h
new file mode 100644
index 0000000..a3dff50
--- /dev/null
+++ b/Imagine/Features/Feat.h
@@ -0,0 +1,54 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+namespace Imagine {
+ /// \addtogroup Features
+ /// @{
+
+ /// alias for feature point coordinates.
+ typedef FVector<float,2> fPoint2;
+
+ /// Feature point.
+ /// Feature point.
+ /// \param T FP descriptor type
+ template <typename T>
+ class FeaturePoint {
+ public:
+ /// FP descriptor type
+ typedef T Desc;
+ /// FP position
+ fPoint2 pos;
+ /// FP scale
+ float scale;
+ /// FP angle
+ float angle;
+ /// FP descriptor
+ Desc desc;
+ /// Read access: FP x position
+ inline float x() const { return pos.x(); }
+ /// Read access: FP y position
+ inline float y() const { return pos.y(); }
+ };
+
+ /// Feature detector.
+ /// Feature detector.
+ /// \param T FeaturePoint type
+ template <typename T>
+ class FeatureDetector {
+ public:
+ /// FP type
+ typedef T Feature;
+ /// Run detection.
+ /// Runs FP detection.
+ /// \param I image on which detection is performed
+ ///
+ /// \dontinclude Features/test/test.cpp \skip detection(
+ /// \skipline runs detection
+ virtual Array<T> run(const Image<byte>& I) const=0;
+ };
+
+ ///@}
+}
diff --git a/Imagine/Features/IO.h b/Imagine/Features/IO.h
new file mode 100644
index 0000000..3991c1f
--- /dev/null
+++ b/Imagine/Features/IO.h
@@ -0,0 +1,209 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+namespace Imagine {
+ /// \addtogroup Features
+ /// @{
+
+ /// Draw on FP.
+ /// Draws a feature point.
+ /// \tparam T FP type
+ /// \param fp FP to draw
+ /// \param off display offset
+ /// \param c color
+ /// \param drawAngle true=draw a radius showing direction
+ /// \param fillIt true=draw a filled circle
+ /// \param z scaling factor (of the whole image)
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline draw one feature
+ template <typename T>
+ void DrawFeature(
+ const T& fp,
+ iPoint2 off=iPoint2(0,0),
+ Color c=Red,
+ bool drawAngle=true,
+ bool fillIt=false,
+ float z=1.f) {
+ if (fillIt)
+ FillCircle(iPoint2(fp.pos*z)+off,int(fp.scale*z+1),c);
+ else
+ DrawCircle(iPoint2(fp.pos*z)+off,int(fp.scale*z+1),c);
+ if (drawAngle)
+ DrawLine(iPoint2(fp.pos*z)+off,
+ iPoint2(fp.pos*z+FVector<float,2>(float(fp.scale*cos(fp.angle)*z),float(fp.scale*sin(fp.angle)*z)))+off,
+ c,1);
+ }
+ /// Draw FPs.
+ /// Draws features points.
+ /// \tparam T FP type
+ /// \param fps array of FPs
+ /// \param off display offset
+ /// \param c color
+ /// \param drawAngle true=draw a radius showing direction
+ /// \param z scaling factor (of the whole image)
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline draw features
+ /// \until ...
+ template <typename T>
+ void DrawFeatures(const Array<T>& fps,
+ iPoint2 off=iPoint2(0,0),
+ Color c=Red,
+ bool drawAngle=true,
+ float z=1.f) {
+
+ NoRefreshBegin();
+ for (size_t i=0;i<fps.size();i++)
+ DrawFeature(fps[i],off,c,drawAngle,false,z);
+ NoRefreshEnd();
+ }
+
+ /// Draw Matches.
+ /// Draws matches between two images, one different (random) color each. Possibility to link matches with a line.
+ /// \tparam T FP type
+ /// \param F1,F2 arrays of FPs
+ /// \param L list of matches
+ /// \param off1,off2 display offsets
+ /// \param z scaling factor (of the whole image)
+ /// \param drawLines true=link matches with lines
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline draw features
+ /// \until ...
+ template <typename T>
+ void DrawMatches(const Array<T>& F1,const Array<T>& F2,
+ const MatchList& L,
+ iPoint2 off1, iPoint2 off2,
+ float z=1.f, bool drawLines=false) {
+ for (MatchList::const_iterator it=L.begin();it!=L.end();++it) {
+ Color c(intRandom(0,255),intRandom(0,255),intRandom(0,255));
+ size_t i1=it->first;
+ size_t i2=it->second;
+ DrawFeature(F1[i1],off1,c,false,true,z);
+ DrawFeature(F2[i2],off2,c,false,true,z);
+ if(drawLines)
+ DrawLine(iPoint2(F1[i1].pos*z)+off1,
+ iPoint2(F2[i2].pos*z)+off2,c,1);
+ }
+ }
+ /// Draw Good Matches.
+ /// Draws matches between two images, good ones in green, bad ones in red. good/bad refers to some given criterion!
+ /// \tparam T FP type
+ /// \tparam C criterion type
+ /// \param F1,F2 arrays of FPs
+ /// \param L list of matches
+ /// \param crit good/bad criterion. \p crit(F1,i1,F2,i2) should be true is OK, false is not
+ /// \param off1,off2 display offsets
+ /// \param z scaling factor (of the whole image)
+ /// \param drawLines true=link matches with lines
+ ///
+ /// \dontinclude Features/test/test.cpp
+ /// \skipline Good/bad criterion based on fundamental matrix
+ /// \until ...
+ /// \skipline draw good matches
+ template <typename T,typename C>
+ void DrawGoodMatches(const Array<T>& F1,const Array<T>& F2,
+ const MatchList& L, C crit,
+ iPoint2 off1, iPoint2 off2,
+ float z=1.f, bool drawLines=false) {
+ for (MatchList::const_iterator it=L.begin();it!=L.end();++it) {
+ size_t i1=it->first;
+ size_t i2=it->second;
+ Color c=crit(F1,i1,F2,i2)?Green:Red;
+ DrawFeature(F1[i1],off1,c,false,true,z);
+ DrawFeature(F2[i2],off2,c,false,true,z);
+ if(drawLines)
+ DrawLine(iPoint2(F1[i1].pos*z)+off1,
+ iPoint2(F2[i2].pos*z)+off2,c,1);
+ }
+ }
+ /// Read FPs.
+ /// Reads feature points to a file.
+ /// \tparam T FP type
+ /// \param feats array of FPs
+ /// \param name file name
+ /// \param loweCompatibility true=Lowe's format (x/y swapped and rounded values)
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline read FP
+ /// \until ...
+ template <typename T>
+ bool ReadFeaturePoints(Array<T>& feats,std::string name,bool loweCompatibility=false) {
+ std::ifstream f(name.c_str());
+ if (!f.is_open()) {
+ std::cerr << "Cant open file" << std::endl;
+ return false;
+ }
+ int nb,sz;
+ f >> nb >> sz;
+ typename T::Desc foo;
+ if (sz!=foo.size()) {
+ std::cerr << "Bad desc size" << std::endl;
+ return false;
+ }
+ feats.resize(nb);
+ for (int i=0;i<nb;i++) {
+ if (loweCompatibility)
+ f >> feats[i].pos[1] >> feats[i].pos[0] >> feats[i].scale >> feats[i].angle;
+ else
+ f >> feats[i].pos[0] >> feats[i].pos[1] >> feats[i].scale >> feats[i].angle;
+ for (int k=0;k<sz;k++) {
+ double x;
+ f >> x;
+ feats[i].desc[k]=typename T::Desc::value_type(x);
+ }
+ }
+ return true;
+ }
+
+
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
+ inline float make_round(float x, int n)
+ {
+ float p=pow(10.f,n);
+ x*=(p);
+ x=(x>=0)?floor(x):ceil(x);
+ return x/p;
+ }
+#endif
+
+ /// Write FPs.
+ /// Writes feature points to a file.
+ /// \tparam T FP type
+ /// \param feats array of FPs
+ /// \param name file name
+ /// \param loweCompatibility true=Lowe's format (x/y swapped and rounded values)
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline write FP
+ /// \until ...
+ template <typename T>
+ bool WriteFeaturePoints(Array<T>& feats,std::string name,bool loweCompatibility=false) {
+ std::ofstream f(name.c_str());
+ if (!f.is_open()) {
+ std::cerr << "Cant open file" << std::endl;
+ return false;
+ }
+ typename T::Desc foo;
+ int sz=foo.size();
+ f << feats.size() << " " << sz << std::endl;
+ for (size_t i=0;i<feats.size();i++) {
+ if (loweCompatibility)
+ f <<make_round(feats[i].y(),2) << ' ' << make_round(feats[i].x(),2) << ' ' << make_round(feats[i].scale,2) << ' ' << make_round(feats[i].angle,3) << std::endl;
+ else
+ f << feats[i].pos[0] << ' ' << feats[i].pos[1] << ' ' << feats[i].scale << ' ' << feats[i].angle << std::endl;
+ for (int k=0;k<sz;k++) {
+ f << double(feats[i].desc[k]) << ' ';
+ if ((k+1)%20==0 || k==sz-1)
+ f << std::endl;
+ }
+ }
+ return true;
+ }
+ ///@}
+
+}
diff --git a/Imagine/Features/Match.h b/Imagine/Features/Match.h
new file mode 100644
index 0000000..51734e6
--- /dev/null
+++ b/Imagine/Features/Match.h
@@ -0,0 +1,106 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+namespace Imagine {
+ /// \addtogroup Features
+ /// @{
+
+ /// Two best matches.
+ /// Find the two best matches for a given FP. Note: distance to second best match is available. Not its index.
+ /// \tparam T FP type
+ /// \param F1,F2 FP arrays to match
+ /// \param i1 index of FP of F1 to match
+ /// \param bi2 index of best match in F2
+ /// \param bd distance between F1[i1].desc and F2[bi2].desc
+ /// \param sbd distance between F1[i1].desc and descriptor of second best match
+ /// \param F pointer to fundamental matrix (optional)
+ /// \param epiDist max distance to epipolar of F1[i1].pos
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline two best matches
+ /// \until ...
+ template <typename T> void find2best(const Array<T> & F1, const Array<T> & F2,
+ size_t i1, size_t & bi2,
+ double & bd, double & sbd,
+ FMatrix<double, 3, 3>* F = 0,
+ double epiDist = 0) {
+ FVector<double,3> m1(F1[i1].pos[0], F1[i1].pos[1], 1);
+ FVector<double,3> l;
+ if (F)
+ l=(*F)*m1;
+
+ bi2 = -1;
+ bd = sbd = std::numeric_limits<double>::max();
+ for (size_t i2 = 0; i2 < F2.size(); i2++) {
+ if (F) {
+ FVector<double,3> m2(F2[i2].pos[0], F2[i2].pos[1], 1);
+ if (abs(l*m2) / std::sqrt(l[0]*l[0] + l[1]*l[1]) > epiDist)
+ continue;
+ }
+
+ double d = distance(F1[i1].desc, F2[i2].desc);
+ if (d < bd) {
+ sbd = bd;
+ bd = d;
+ bi2 = i2;
+ } else if (d < sbd) {
+ sbd = d;
+ }
+ }
+
+ }
+
+ /// Alias to list of pairs of indices
+ typedef std::list<std::pair<size_t,size_t> > MatchList;
+
+ /// Lowe's match.
+ /// Matching based on Lowe's recommendations: ratio between best and second best matches. Naive impl. without ANN
+ /// \tparam T FP type
+ /// \param F1,F2 FP arrays to match
+ /// \param ratio First/Second ratio
+ /// \param sym symmetric test: match OK if OK both ways
+ /// \param descDist maximum descriptor distance (0=none)
+ /// \param F pointer to fundamental matrix (optional)
+ /// \param epiDist max distance to epipolar of F1[i1].pos
+ /// \return list of matches
+ ///
+ /// \dontinclude Features/test/test.cpp \skip matching()
+ /// \skipline Lowe's match
+ /// \until ...
+ template <typename T> MatchList LoweMatch(const Array<T>& F1,
+ const Array<T>& F2,
+ double ratio = 0.6,
+ bool sym = false,
+ double descDist = 0,
+ FMatrix<double, 3, 3>* F = 0,
+ double epiDist = 0) {
+ MatchList L;
+ FMatrix<double, 3, 3> FTrans; // F transpose
+ FMatrix<double, 3, 3>* FT=0; // pointer (or 0)
+ if (F) {
+ FTrans = transpose(*F);
+ FT = &FTrans;
+ }
+ for (size_t i1=0;i1<F1.size();i1++) {
+ size_t bi2;
+ double bd,sbd;
+ find2best(F1,F2,i1,bi2,bd,sbd,F,epiDist);
+ if (bi2!=-1 && bd/sbd < ratio && (descDist==0 || bd <descDist)) {
+ if (!sym)
+ L.push_back(std::pair<size_t,size_t>(i1,bi2));
+ else {
+ size_t bi1;
+ find2best(F2,F1,bi2,bi1,bd,sbd,FT,epiDist);
+ if (bd/sbd<ratio && bi1==i1) // i1 -> bi2 -> bi1 with bi1==i1
+ L.push_back(std::pair<size_t,size_t>(i1,bi2));
+ }
+ }
+ }
+ return L;
+ }
+ ///@}
+
+}
diff --git a/Imagine/Features/SIFT.h b/Imagine/Features/SIFT.h
new file mode 100644
index 0000000..4f39fca
--- /dev/null
+++ b/Imagine/Features/SIFT.h
@@ -0,0 +1,77 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+// SIFT from VLFeat
+
+namespace Imagine {
+ /// \addtogroup Features
+ /// @{
+
+ /// SIFT descriptor
+ typedef FeaturePoint<FVector<byte,128> > SIFT;
+
+ /// SIFT detector. VLFeat implementation.
+ /// SIFT detector. VLFeat implementation.
+ class SIFTDetector:public FeatureDetector<SIFT> {
+ // First Octave Index.
+ int firstOctave;
+ // Number of octaves.
+ int numOctaves;
+ // Number of scales per octave.
+ int numScales;
+ // Max ratio of Hessian eigenvalues.
+ float edgeThresh;
+ // Min contrast.
+ float peakThresh;
+ public:
+ /// Constructor.
+ /// Constructor.
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, constructor
+ SIFTDetector() {
+ firstOctave=-1;
+ numOctaves=-1;
+ numScales=3;
+ edgeThresh=10.0f; peakThresh=0.04f;
+ }
+ /// Number of octaves.
+ /// Sets number of octaves. -1 = max. default=-1
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, number of octaves
+ void setNumOctaves(int n) { numOctaves=n; }
+ /// First octave.
+ /// Sets first octave. -1 = doubleImageSize. default=-1
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, first octave
+ void setFirstOctave(int n) { firstOctave=n; }
+ /// Number of scales per octave.
+ /// Sets number of scales per octave. default=3
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, number of scales per octave
+ void setNumScales(int n) { numScales=n; }
+ /// Max ratio of Hessian eigenvalues.
+ /// Sets max ratio of Hessian eigenvalues. default=10
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, max ratio of Hessian eigenvalues
+ void setEdgeThresh(float t) { edgeThresh=t; }
+ /// Min contrast.
+ /// Sets min contrast. Scale between 0 and 1. Default=0.04
+ ///
+ /// \dontinclude Features/test/test.cpp \skip main()
+ /// \skipline sift, min contrast
+ void setPeakThresh(float t) { peakThresh=t; }
+
+ // Implementation
+ Array<SIFT> run(const Image<byte>& I) const;
+ };
+ ///@}
+
+}
diff --git a/Imagine/Optim.h b/Imagine/Optim.h
new file mode 100644
index 0000000..5e5d7c3
--- /dev/null
+++ b/Imagine/Optim.h
@@ -0,0 +1,20 @@
+// ===========================================================================
+// Imagine++ Libraries
+// Copyright (C) Imagine
+// For detailed information: http://imagine.enpc.fr/software
+// ===========================================================================
+
+#ifndef _IMAGINEOPTIM_H
+#define _IMAGINEOPTIM_H
+
+#include <Imagine/Common.h>
+
+/// \defgroup Optim Optim Library
+/// @{
+/// \example Optim/test/example.cpp
+/// @}
+
+#include "Optim/Ransac.h" // Ransac
+
+#endif
+
diff --git a/Imagine/Optim/Ransac.h b/Imagine/Optim/Ransac.h
new file mode 100644
index 0000000..2ad56c7
--- /dev/null
+++ b/Imagine/Optim/Ransac.h
@@ -0,0 +1,207 @@
+/// \file
+/// \brief Robust estimation in presence of outliers with RANSAC and variants.
+/// \details Two functions are proposed: Ransac (fix threshold for outliers and
+/// minimize number of outliers) and LMedS (define outliers as quantile of worst
+/// residuals and minimize the threshold). Template parameters should statisfy:
+/// \li \c SampleSize number of data samples necessary to evaluate parameters.
+/// \li \c DataInputIterator a forward iterator type over data.
+/// \li \c Estimator functor of arguments
+/// (FArray<DataInputIterator,SampleSize>,OutputIterator) putting in
+/// \c OutputIterator the sets of parameters evaluated from samples in array
+/// argument. OutputIterator::value_type must be \c Parameters.
+/// \li \c Residual functor of arguments
+/// (Parameters,DataInputIterator::value_type) returning a \c double error of
+/// sample relative to the parameters.
+/// \li \c Parameters set of parameters to estimate.
+///
+/// Normally, only the \c SampleSize template parameter must be precised, the
+/// others are implicitely deduced from the arguments by the compiler.
+#include <vector>
+
+namespace Imagine {
+
+ /// \addtogroup Optim
+ /// @{
+
+ /// \brief Get a random set of \a SampleSize samples out of \a n elements.
+ /// \details \a out iterates over collection of type \a InputIterator.
+ /// Only one pass is done through \a first, so it remains efficient in
+ /// case of a forward iterator. If \a bUnique is set and the samples are not
+ /// unique, return \a false. This is unlikely if SampleSize<<n, but in case the
+ /// best is to try again.
+ inline size_t size_tRandom(size_t a)
+ {
+ return (size_t)((a+.999)*doubleRandom());
+ }
+ template <int SampleSize, class InputIterator, class OutputIterator>
+ bool getSample(InputIterator first, size_t n, OutputIterator out,
+ bool bUnique=false)
+ {
+ FArray<std::pair<size_t,int>,SampleSize> index;
+ for(int i=0; i < SampleSize; i++)
+ index[i] = std::make_pair(size_tRandom(n-1), i);
+ std::sort(index.begin(), index.end()); // Lexicographic order
+ size_t prev=0;
+ for(int i=0; i < SampleSize; i++) {
+ if(bUnique && i != 0 && index[i].first == prev)
+ return false;
+ std::advance(first, index[i].first-prev);
+ prev = index[i].first;
+ OutputIterator o=out;
+ std::advance(o, index[i].second);
+ *o = first;
+ }
+ return true;
+ }
+
+ /// Number of samplings to have at least \a minProba probability of absence of
+ /// outlier in a sample of \a SampleSize elements.
+ inline int getNumSamples(double minProba, double outlierRatio, int SampleSize)
+ {
+ return int( std::log(1.-minProba) /
+ std::log(1.-std::pow(1.-outlierRatio, SampleSize)) );
+ }
+
+ /// \brief The famous Random Sample Consensus algorithm (Fischler&Bolles 1981).
+ /// \details The number of tests is reevaluated down as soon as more inliers are
+ /// found. Return the number of inliers.
+ template <int SampleSize, class DataInputIterator, class Estimator,
+ class Residual, class Parameters>
+ std::size_t Ransac(DataInputIterator first, DataInputIterator last,
+ const Estimator& estimator, const Residual& residual,
+ Parameters& bestParameters,
+ double outlierThreshold,
+ double maxOutlierRatio=0.5,
+ double minProba=0.99)
+ {
+ std::size_t bestInliers = 0;
+ std::size_t n = std::distance(first,last);
+
+ FArray<DataInputIterator,SampleSize> sample;
+
+ // Required number of iterations is evaluated from outliers ratio
+ int N = (SampleSize<=n)?
+ getNumSamples(minProba, maxOutlierRatio, SampleSize): 0;
+ for (int i=0; i < N; i++) {
+ // Try until getting a set of SampleSize unique samples
+ while(! getSample<SampleSize>(first, n, sample.begin(), true)) {}
+ // Estimate parameters: the solutions are stored in a vector
+ std::vector<Parameters> paramVec;
+ estimator(sample, std::back_inserter(paramVec));
+
+ // Now test the solutions on the whole data
+ typename std::vector<Parameters>::iterator it = paramVec.begin();
+ for(; it != paramVec.end(); ++it) {
+ std::size_t inliers = 0;
+ for(DataInputIterator it2=first; it2 != last; ++it2)
+ if(residual(*it,*it2) < outlierThreshold)
+ ++ inliers;
+
+ // Store best solution and reestimate number of tests
+ if (inliers > bestInliers) {
+ bestInliers = inliers;
+ bestParameters = *it;
+ double outlierRatio = (n - inliers)/(double)n;
+ if(outlierRatio < maxOutlierRatio)
+ N = getNumSamples(minProba, outlierRatio, SampleSize);
+ }
+ }
+ }
+ return bestInliers;
+ }
+
+ /// \brief Variant of RANSAC using Least Median of Squares.
+ /// \details Instead of using a fixed threshold to distinguish inlier/outlier,
+ /// find the threshold at 1-\a outlierRatio quantile of residuals and keep the
+ /// parameters minimizing this threshold. The final threshold
+ /// returned in \a outlierThreshold is a multiple of this and can be used to
+ /// filter out outliers.
+ template <int SampleSize, class DataInputIterator,
+ class Estimator, class Residual, class Parameters>
+ double LMedS(DataInputIterator first, DataInputIterator last,
+ const Estimator& estimator, const Residual& residual,
+ Parameters& bestParameters,
+ double* outlierThreshold=NULL,
+ double outlierRatio=0.5,
+ double minProba=0.99)
+ {
+ std::size_t n = std::distance(first,last); // Data size
+ std::vector<double> residuals(n); // Array for storing residuals
+
+ FArray<DataInputIterator,SampleSize> sample;
+ double bestMedian = std::numeric_limits<double>::max();
+
+ // Required number of iterations is evaluated from outliers ratio
+ const int N = (SampleSize<n)?
+ getNumSamples(minProba, outlierRatio, SampleSize): 0;
+ for (size_t i=0; i < N; i++) {
+ // Try until getting a set of SampleSize unique samples
+ while(! getSample<SampleSize>(first, n, sample.begin(), true)) {}
+ // Estimate parameters: the solutions are stored in a vector
+ std::vector<Parameters> paramVec;
+ estimator(sample, std::back_inserter(paramVec));
+
+ // Now test the solutions on the whole data
+ typename std::vector<Parameters>::iterator it = paramVec.begin();
+ for(; it != paramVec.end(); ++it) {
+ // Compute residuals
+ std::vector<double>::iterator out = residuals.begin();
+ for(DataInputIterator it2 = first; it2 != last; ++it2)
+ *out++ = residual(*it,*it2);
+
+ // Compute median
+ std::vector<double>::iterator itMedian = residuals.begin() +
+ std::size_t( n*(1.-outlierRatio) );
+ std::nth_element(residuals.begin(), itMedian, residuals.end());
+ double median = *itMedian;
+
+ // Store best solution
+ if(median < bestMedian) {
+ bestMedian = median;
+ bestParameters = *it;
+ }
+ }
+ }
+
+ /* this array of precomputed values corresponds to the inverse
+ cumulative function for a normal distribution. For more information
+ consult the litterature (Robust Regression for Outlier Detection,
+ rouseeuw-leroy). The values are computed for each 5% */
+ static const double ICDF[21] = {
+ 1.4e16, 15.94723940, 7.957896558, 5.287692054,
+ 3.947153876, 3.138344200, 2.595242369, 2.203797543,
+ 1.906939402, 1.672911853, 1.482602218, 1.323775627,
+ 1.188182950, 1.069988721, 0.9648473415, 0.8693011162,
+ 0.7803041458, 0.6946704675, 0.6079568319,0.5102134568,
+ 0.3236002672
+ };
+
+ // Evaluate the outlier threshold
+ if(outlierThreshold) {
+ double sigma = ICDF[int((1.-outlierRatio)*20.)] *
+ (1. + 5. / double(n - SampleSize));
+ *outlierThreshold = (double)(sigma * sigma * bestMedian * 4.);
+ }
+
+ return bestMedian;
+ }
+
+ /// \brief Keep only inliers in data set \a first to \a last.
+ /// \details The returned iterator marks the new end of the list. For example,
+ /// if data is in a \code std::list ls \endcode, you can call
+ /// \code ls.erase(FilterOutliers(ls.begin(),ls.end(),res,param,th), ls.end());
+ /// \endcode
+ template <class DataIterator, class Residual, class Parameters>
+ DataIterator FilterOutliers(DataIterator first, DataIterator last,
+ Residual residual,
+ Parameters param, double outlierThreshold)
+ {
+ DataIterator it=first;
+ for(; first != last; first++)
+ if(residual(param,*first) < outlierThreshold)
+ *it++ = *first;
+ return it;
+ }
+ ///@}
+
+} \ No newline at end of file
diff --git a/display.cpp b/display.cpp
new file mode 100644
index 0000000..62efefc
--- /dev/null
+++ b/display.cpp
@@ -0,0 +1,37 @@
+#include "display.h"
+
+ImageSIFT MergePair(const ImageSIFT & I1, const ImageSIFT & I2, const Match & M) {
+ int w = I1.width() + I2.width(), h = I1.height() + I2.height();
+ Image<Color> Merged(w, h);
+ Merged.fill(Black);
+ int i, j;
+ double n;
+ Vec2 t;
+ RGB<double> colorTemp;
+ for (j = - h/4; j < 3*h/4; j++) {
+ for (i = - w/4; i < 3*w/4; i++) {
+ n = 0.;
+ colorTemp = RGB<double>(0., 0., 0.);
+ t = applyHomography(M.H, Vec2(i, j));
+ if (i >= 0 && j >= 0 && j < I1.height() && i < I1.width()) {
+ colorTemp = RGB<double>(I1(i, j));
+ n += 1.;
+ }
+ if (t[0] >= 0 && t[0] < I2.width() && t[1] >= 0 && t[1] < I2.height()) {
+ colorTemp += RGB<double>(I2(t[0], t[1]));
+ n += 1.;
+ }
+ if (n > 0)
+ Merged(i + w/4, j + h/4) = Color((1. / n) * colorTemp);
+ else
+ Merged(i + w/4, j + h/4) = Color(0, 0, 0);
+ }
+ }
+ return ImageSIFT(Merged);
+}
+
+void DisplayImage(const ImageSIFT & I) {
+ Window W = OpenWindow(I.width(), I.height());
+ SetActiveWindow(W);
+ display(I.get_data());
+}
diff --git a/display.h b/display.h
new file mode 100644
index 0000000..90b346e
--- /dev/null
+++ b/display.h
@@ -0,0 +1,14 @@
+#ifndef __DISPLAY_H_191109__
+#define __DISPLAY_H_191109__
+
+#include "matches.h"
+#include "images.h"
+#include <Imagine/Images.h>
+using namespace std;
+using namespace Imagine;
+
+ImageSIFT MergePair(const ImageSIFT & I1, const ImageSIFT & I2,
+ const Match & M);
+void DisplayImage(const ImageSIFT & I);
+
+#endif
diff --git a/image0006.jpg b/image0006.jpg
new file mode 100644
index 0000000..293f0ea
--- /dev/null
+++ b/image0006.jpg
Binary files differ
diff --git a/image0007.jpg b/image0007.jpg
new file mode 100644
index 0000000..b45833a
--- /dev/null
+++ b/image0007.jpg
Binary files differ
diff --git a/image0008.jpg b/image0008.jpg
new file mode 100644
index 0000000..78bfc5d
--- /dev/null
+++ b/image0008.jpg
Binary files differ
diff --git a/images.cpp b/images.cpp
new file mode 100644
index 0000000..329e5ae
--- /dev/null
+++ b/images.cpp
@@ -0,0 +1,17 @@
+#include "images.h"
+#include "framework.h"
+
+ImageSIFT::ImageSIFT (string path) {
+ load(this->data, strSRCPATH(path));
+ this->ComputeSifts();
+}
+
+ImageSIFT::ImageSIFT (const Image<Color> & source)
+ :data(source) {
+ this->ComputeSifts();
+}
+
+void ImageSIFT::ComputeSifts () {
+ Detector d;
+ this->sifts = d.run(this->data);
+}
diff --git a/images.h b/images.h
new file mode 100644
index 0000000..20d6b19
--- /dev/null
+++ b/images.h
@@ -0,0 +1,33 @@
+#ifndef __IMAGES_H_191109__
+#define __IMAGES_H_191109__
+
+#include <string>
+using namespace std;
+
+#include <Imagine/Images.h>
+#include "Imagine/Features.h"
+using namespace Imagine;
+
+//TODO un truc moins immonde
+#ifndef FEAT_DE
+#define FEAT_DE
+typedef SIFTDetector Detector;
+typedef Detector::Feature Feat;
+#endif
+
+class ImageSIFT {
+private:
+ Image<Color> data;
+ Array<Feat> sifts;
+ void ComputeSifts ();
+public:
+ ImageSIFT (string path);
+ ImageSIFT (const Image<Color> & source);
+ inline const Array<Feat> & getSIFTS() const {return this->sifts;};
+ inline int width() const {return this->data.width();};
+ inline int height() const {return this->data.height();};
+ inline Color operator () (int i, int j) const {return this->data(i, j);};
+ inline const Image<Color> & get_data () const {return this->data;};
+};
+
+#endif
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..b0dde00
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,54 @@
+#include <Imagine/Images.h>
+#include <Imagine/LinAlg.h>
+
+#include "matches.h"
+#include "images.h"
+#include "display.h"
+#include <vector>
+#include <string>
+#include <iostream>
+
+using namespace std;
+using namespace Imagine;
+
+int main() {
+ /*
+#if 0
+string imnames[] = {"image0006.jpg", "image0007.jpg", "image0008.jpg"};
+ const int n_images = 3;
+#else
+ string imnames[] = {"image0006.jpg", "image0007.jpg"};
+ const int n_images = 2;
+#endif
+ vector<ImageSIFT> images;
+ for (int i = 0; i < n_images; i++)
+ images.push_back(ImageSIFT(imnames[i]));
+
+ vector<Match> m = ComputeAllMatches(images);
+ for (int i = 0; i < m.size(); i++) {
+ cout << "matches : " << m[i].matches
+ << " inliers : " << m[i].inliers << endl;
+ }
+ DisplayPair(images[0], images[1], m[0]);
+ Terminate();
+ */
+ string imnames[] = {"image0006.jpg", "image0007.jpg", "image0008.jpg"};
+ const int n_images = 3;
+ vector<ImageSIFT> images;
+ for (int i = 0; i < n_images; i++)
+ images.push_back(ImageSIFT(imnames[i]));
+
+ Match m12 = ComputeMatches(images[0], images[1]);
+ cout << "12 computed" << endl;
+ ImageSIFT i12 = MergePair(images[1], images[1], m12);
+ cout << "12 merged" << endl;
+ Match m123 = ComputeMatches(i12, images[2]);
+ cout << "123 computed" << endl;
+ ImageSIFT i123 = MergePair(i12, images[2], m123);
+ cout << "123 merged" << endl;
+ DisplayImage(i123);
+ Click();
+ Terminate();
+
+}
+
diff --git a/matches.cpp b/matches.cpp
new file mode 100644
index 0000000..8b1429d
--- /dev/null
+++ b/matches.cpp
@@ -0,0 +1,28 @@
+#include "matches.h"
+#include <algorithm>
+
+Match ComputeMatches (const ImageSIFT & I1, const ImageSIFT & I2) {
+ Match ret;
+ MatchList matches = LoweMatch(I1.getSIFTS(),I2.getSIFTS(), 0.5, true);
+ ret.matches = matches.size();
+ cout << matches.size() << endl;
+ Ransac(I1.getSIFTS(), I2.getSIFTS(), matches);
+ ret.inliers = matches.size();
+ cout << matches.size() << endl;
+ ret.H = meanHomography(I1.getSIFTS(), I2.getSIFTS(), matches);
+ return ret;
+}
+
+vector<Match> ComputeAllMatches (const vector<ImageSIFT> & images) {
+ //TODO : visiblement buggé
+ vector<Match> ret;
+ int n = images.size();
+ int i, j;
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ if (i != j) {
+ ret.push_back(ComputeMatches(images[i], images[j]));
+ }
+ return ret;
+ //sort(ret.begin(), ret.end());
+}
diff --git a/matches.h b/matches.h
new file mode 100644
index 0000000..649c507
--- /dev/null
+++ b/matches.h
@@ -0,0 +1,21 @@
+#ifndef __MATCHES_H_191109__
+#define __MATCHES_H_191109__
+
+#include <vector>
+using namespace std;
+
+#include "framework.h"
+#include "images.h"
+
+struct Match {
+ int matches, inliers;
+ Homography H;
+ inline bool operator < (const Match & source) const {
+ return this->inliers < source.inliers;
+ }
+};
+
+Match ComputeMatches (const ImageSIFT & I1, const ImageSIFT & I2);
+vector<Match> ComputeAllMatches (const vector<ImageSIFT> & v);
+
+#endif
diff --git a/vl/generic.c b/vl/generic.c
new file mode 100644
index 0000000..48de409
--- /dev/null
+++ b/vl/generic.c
@@ -0,0 +1,400 @@
+/** @file generic.c
+ ** @author Andrea Vedaldi
+ ** @brief Generic - Definition
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+/**
+ @mainpage VLFeat -- Vision Lab Features Library
+
+ @version __VLFEAT_VERSION__
+ @author Andrea Vedaldi (vedaldi@cs.ucla.edu)
+ @author Brian Fulkerson (bfulkers@cs.ucla.edu)
+
+ @par Copyright &copy; 2007-09 Andrea Vedaldi and Brian Fulkerson
+
+ <em>VLFeat C library contains implementations of common computer
+ vision algorithms, with a special focus on visual features for
+ matching image regions. Applications include structure from motion
+ and object and category detection and recognition.
+
+ We strive to make the library free of clutter, portable (VLFeat is
+ largely C-89 compatible), and self- documented. Different parts of
+ the library are weakly interdependent, simplifying understanding and
+ extraction of code.</em>
+
+ @section main-contents Contents
+
+ - @ref design "Design Concepts"
+ - @ref design-objects "Objects"
+ - @ref design-resources "Memory and resource management"
+ - @ref design-threads "Multi-threading"
+ - @ref design-portability "Portability"
+ - @ref main-glossary "Glossary"
+ - Support functionalities
+ - @ref host.h "Platform abstraction"
+ - @ref generic.h "Errors, memory, logging, and others"
+ - @ref mathop.h "Math operations"
+ - @ref random.h "Random number generator"
+ - @ref stringop.h "String operations"
+ - @ref imop.h "Image operations"
+ - @ref rodrigues.h "Rodrigues formula"
+ - @ref mexutils.h "MATLAB MEX helper functions"
+ - Algorithms
+ - @ref sift.h "Scale Invariant Feature Transform (SIFT)"
+ - @ref dsift.h "Dense SIFT (DSIFT)"
+ - @ref mser.h "Maximally Stable Extremal Regions (MSER)"
+ - @ref ikmeans.h "Integer K-means (IKM)"
+ - @ref hikmeans.h "Hierarchical Integer K-means (HIKM)"
+ - @ref aib.h "Agglomerative Information Bottleneck (AIB)"
+ - @ref quickshift.h "Quick shift image segmentation"
+
+ @section design VLFeat Design Concepts
+
+ VLFeat follows a simple but rigorous design that makes it portable
+ and simple to use in conjunction with high level language such as
+ MATLAB. This section illustrates and motivates the aspects of the
+ design that are relevant to the users of the library. Most of the
+ features discussed in here are implemented in the @ref generic.h
+ module.
+
+ @subsection design-objects Objects
+
+ Most of VLFeat functionalities are implemented as opaque data
+ structures, to which we refer as "objects". Typically, you create an
+ object by means of a constructor function and dispose of it by means
+ of a destructor function. The data member of the object should not be
+ accessed directly, but by means of appropriate accessor functions
+ (typically containing the @c _get and _set keywords in their names).
+
+ @subsection design-resources Memory and Resource Management
+
+ Resource management in VLFeat is minimal. In most cases, you can
+ assume that VLFeat does not provide any resource management
+ functionality at all. Objects or memory blocks allocated by the
+ library but owned by the client must be explicitly disposed. The
+ following rule helps identifying such blocks and objects:
+
+ <b> The client owns a memory block or object if, and only if, it is
+ returned by a library call containing the keywords @c _new or @c
+ _copy, or by the allocator functions ::vl_malloc, ::vl_calloc,
+ ::vl_realloc.</b>
+
+ More in detail, the following rules apply:
+
+ - Memory is the only managed resource. Other resources used by the
+ library (e.g. files) are either passed to the library by the
+ client, or acquired and released within the scope of a library
+ function call.
+ - The memory manager can be customized through ::vl_set_alloc_func
+ (which sets the implementations of ::vl_malloc, ::vl_realloc,
+ ::vl_calloc and ::vl_free). The library allocates memory only through
+ these functions.
+ - The memory manager is global to all threads.
+ - At any moment, there is only one memory manager in
+ existence. ::vl_set_alloc_func can be used only before any other
+ function is invoked, or right before or after the previous memory
+ manager has relinquished all memory.
+ - Such rules apply both to the library client and to the library
+ implementation. The only exception regards the memory allocated in
+ the global library state, which uses the native memory manager.
+
+ These rules make it possible to undo all the work done by the
+ library at any given instant. Disposing the memory allocated by the
+ custom memory manager essentially "resets" the library (except for
+ its global state). This functionality is used extensively in the
+ implementation of MATLAB MEX functions that access the library to
+ support abrupt interruption.
+
+ @note The rules guarantee that all memory allocated by the library
+ at any given time is allocated by the same memory manager, except
+ for the global state that should survive the "reset"
+ operation. In order to support multiple memory managers, one
+ should keep track of the allocator of each object (or memory
+ block). Moreover, partial de-allocation of a pool of objects is
+ dangerous, as such objects may be referred by other objects that
+ are not being de-allocated, corrupting their state. A proper
+ solution to the latter problem is the retain/release mechanism
+ implemented, for instance, by Apple Core Foundation or Cocoa.
+
+ @subsection design-threads Multi-threading
+
+ The library is currently <em>not</em> thread safe, but this support
+ will be added in a future release.
+
+ The library is almost entirely reentrant. The only thread-sensitive
+ operations are on the global library state and are limited to:
+
+ - Global library configuration (e.g. ::vl_set_alloc_func).
+ - Random number generator state (@ref random.h).
+ - Error handling (e.g. ::vl_err_no).
+
+ @subsection design-portability Portability features
+
+ Platform dependent details are isolated in the @ref generic.h
+ library module. These include:
+
+ - Atomic types (e.g. ::vl_int32).
+ - Special syntaxes for the declaration of symbols exported by the library
+ and inline functions (e.g. ::VL_EXPORT).
+ - Host-dependent conversion of data endianess
+ (e.g. ::vl_swap_host_big_endianness_8()).
+
+ VLFeat uses processor specific features (e.g. Intel SSE) if those
+ are available at compile time.
+
+ <!-- @see http://www.macresearch.org/how_to_properly_use_sse3_and_ssse3_and_future_intel_vector_extensions_0 -->
+
+ @section main-glossary Glossary
+
+ - <b>Column-major.</b> A <em>M x N </em> matrix <em>A</em> is
+ stacked with column-major order as the sequence \f$(A_{11}, A_{21},
+ \dots, A_{12}, \dots)\f$. More in general, when stacking a multi
+ dimensional array this indicates that the first index is the one
+ varying most quickly, with the other followed in the natural order.
+ - <b>Opaque structure.</b> A structure is opaque if the user is not supposed
+ to access its member directly, but through appropriate interface functions.
+ Opaque structures are commonly used to define objects.
+ - <b>Row-major.</b> A <em>M x N </em> matrix <em>A</em> is
+ stacked with row-major order as the sequence \f$(A_{11}, A_{12},
+ \dots, A_{21}, \dots)\f$. More in general, when stacking a multi
+ dimensional array this indicates that the last index is the one
+ varying most quickly, with the other followed in reverse order.
+ - <b>Feature frame.</b> A <em>feature frame</em> is the geometrical
+ description of a visual features. For instance, the frame of
+ a @ref sift.h "SIFT feature" is oriented disk and the frame of
+ @ref mser.h "MSER feature" is either a compact and connected set or
+ a disk.
+ - <b>Feature descriptor.</b> A <em>feature descriptor</em> is a quantity
+ (usually a vector) which describes compactly the appearance of an
+ image region (usually corresponding to a feature frame).
+**/
+
+/**
+ @file generic.h
+
+ This module provides the following functionalities:
+
+ - @ref generic-preproc
+ - @ref generic-error
+ - @ref generic-heap
+ - @ref generic-logging
+ - @ref generic-time
+
+ @section generic-preproc C preprocessor
+
+ VLFeat provides a few C preprocessor macros of general
+ utility. These include stringification (::VL_STRINGIFY,
+ ::VL_XSTRINGIFY) and concatenation (::VL_CAT, ::VL_XCAT) of symbols.
+
+ @section generic-error Error handling
+
+ Some VLFeat functions signal errors as the standard C library. Such
+ functions return 0 when they succeed and -1 when they fail, and set
+ the global variable ::vl_err_no with a code identifying the error
+ occurred. This variable is never set on success and should be
+ examined right after an error had occurred.
+
+ @section generic-heap Heap allocation
+
+ VLFeat uses the ::vl_malloc(), ::vl_realloc(), ::vl_calloc() and
+ ::vl_free() functions to allocate the heap. Normally these functions
+ are mapped to the underlying standard C library implementations. However
+ ::vl_set_alloc_func() can be used to map them to other implementations.
+ For instance, in MATALB MEX files these functions are mapped to
+ the MATLAB equivalent which has a garbage collection mechanism to cope
+ with interruptions during execution.
+
+ @section generic-logging Logging
+
+ VLFeat uses the macros ::VL_PRINT and ::VL_PRINTF to print progress or
+ debug informations. These functions are normally mapped to the @c
+ printf function of the underlying standard C library. However
+ ::vl_set_printf_func() can be used to map it to a different
+ implementation. For instance, in MATLAB MEX files this function is
+ mapped to @c mexPrintf. Setting the function to @c NULL disables
+ logging.
+
+ @section generic-time Measruing time
+
+ VLFeat provides ::vl_tic() and ::vl_toc() as an easy way of
+ measuring elapsed time.
+
+**/
+
+#include "generic.h"
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#ifdef VL_OS_WIN
+#include <Windows.h>
+#endif
+
+VL_EXPORT int vl_err_no = 0 ;
+VL_EXPORT char vl_err_msg [VL_ERR_MSG_LEN + 1] = "" ;
+
+/** ------------------------------------------------------------------
+ ** @brief Get version string
+ ** @return library version string
+ **/
+
+VL_EXPORT
+char const * vl_get_version_string ()
+{
+ return VL_VERSION_STRING ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Print information about the library
+ ** @return library version string
+ **/
+
+VL_EXPORT
+void vl_print_info ()
+{
+ VL_PRINTF ("VLFeat version %s\n", vl_get_version_string()) ;
+ vl_print_host_info () ;
+}
+
+/** @internal@brief A printf that does not do anything */
+static int
+do_nothing_printf (char const* format, ...)
+{
+ return 0 ;
+}
+
+/** @internal@brief Customizable @c malloc function pointer */
+void *(*vl_malloc_func) (size_t) = &malloc ;
+
+/** @internal@brief Customizable @c realloc function pointer */
+void *(*vl_realloc_func) (void*,size_t) = &realloc ;
+
+/** @internal@brief Customizable @c calloc function pointer */
+void *(*vl_calloc_func) (size_t, size_t) = &calloc ;
+
+/** @internal@brief Customizable @c free function pointer */
+void (*vl_free_func) (void*) = &free ;
+
+/** @internal@brief Customizable @c printf function pointer */
+int (*vl_printf_func) (char const *, ...)= printf ; /* &do_nothing_printf ;*/
+
+/** --------------------------------------------------------------- */
+
+/** @fn ::vl_malloc(size_t)
+ ** @brief Call customizable @c malloc function
+ ** @param n number of bytes to allocate.
+ **
+ ** The function calls the user customizable @c malloc.
+ **
+ ** @return result of @c malloc
+ **/
+
+/** @fn ::vl_realloc(void*,size_t)
+ ** @brief Call customizable @c resize function
+ **
+ ** @param ptr buffer to reallocate.
+ ** @param n number of bytes to allocate.
+ **
+ ** The function calls the user-customizable @c realloc.
+ **
+ ** @return result of the user-customizable @c realloc.
+ **/
+
+/** @fn ::vl_calloc(size_t,size_t)
+ ** @brief Call customizable @c calloc function
+ **
+ ** @param n size of each element in byte.
+ ** @param size size of the array to allocate (number of elements).
+ **
+ ** The function calls the user-customizable @c calloc.
+ **
+ ** @return result of the user-customizable @c calloc.
+ **/
+
+/** @fn ::vl_free(void*)
+ ** @brief Call customizable @c free function
+ **
+ ** @param ptr buffer to free.
+ **
+ ** The function calls the user customizable @c free.
+ **/
+
+/** ------------------------------------------------------------------
+ ** @brief Set memory allocation functions
+ ** @param malloc_func pointer to @c malloc.
+ ** @param realloc_func pointer to @c realloc.
+ ** @param calloc_func pointer to @c calloc.
+ ** @param free_func pointer to @c free.
+ **/
+
+VL_EXPORT
+void vl_set_alloc_func (void *(*malloc_func) (size_t),
+ void *(*realloc_func) (void*, size_t),
+ void *(*calloc_func) (size_t, size_t),
+ void (*free_func) (void*))
+{
+ vl_malloc_func = malloc_func ;
+ vl_realloc_func = realloc_func ;
+ vl_calloc_func = calloc_func ;
+ vl_free_func = free_func ;
+}
+
+VL_EXPORT
+void
+vl_set_printf_func (printf_func_t printf_func)
+{
+ vl_printf_func = printf_func ? printf_func : do_nothing_printf ;
+}
+
+#ifdef VL_OS_WIN
+LARGE_INTEGER tic_freq ;
+LARGE_INTEGER tic_mark ;
+#else
+clock_t tic_mark ; /**< @internal Store clock time for ::vl_tic() */
+#endif
+
+/** ------------------------------------------------------------------
+ ** @brief Set time reference
+ **/
+
+void vl_tic()
+{
+#ifdef VL_OS_WIN
+ QueryPerformanceFrequency(&tic_freq) ;
+ QueryPerformanceCounter(&tic_mark) ;
+#else
+ tic_mark = clock() ;
+#endif
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get time since reference
+ **
+ ** Returns the processor time elapsed since ::vl_tic() was called.
+ **
+ ** @remark On UNIX, this function uses the @c clock() system call.
+ ** On Windows, it uses the @c QueryPerformanceCounter() system call,
+ ** which is more accurate than @c clock() on this platform.
+ **
+ ** @return time in seconds.
+ **/
+
+double vl_toc()
+{
+#ifdef VL_OS_WIN
+ LARGE_INTEGER toc_mark ;
+ QueryPerformanceCounter(&toc_mark) ;
+ return (double) (toc_mark.QuadPart - tic_mark.QuadPart) / tic_freq.QuadPart ;
+#else
+ return (double) (clock() - tic_mark) / CLOCKS_PER_SEC ;
+#endif
+}
diff --git a/vl/generic.h b/vl/generic.h
new file mode 100644
index 0000000..9262c7a
--- /dev/null
+++ b/vl/generic.h
@@ -0,0 +1,255 @@
+/** @file generic.h
+ ** @author Andrea Vedaldi
+ ** @brief Generic
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#ifndef VL_GENERIC_H
+#define VL_GENERIC_H
+
+#include <stddef.h>
+#include <time.h>
+#include <assert.h>
+
+#include "host.h"
+
+/** @brief Library version string */
+#define VL_VERSION_STRING "0.9.4"
+
+/** ------------------------------------------------------------------
+ ** @name C preprocssor
+ ** @{ */
+
+/** @brief Convert the argument to a string
+ **
+ ** @param x value to be stringified.
+ **
+ ** This macro stringifies the argument @a x by means of the
+ ** <code>#</code> prerpocessor operator.
+ **
+ ** The standard C preprocessor does not prescan arguments which are
+ ** stringified, so
+ **
+ ** @code
+ ** #define A B
+ ** char const * str = VL_STRINGIFY(A) ;
+ ** @endcode
+ **
+ ** initializes <code>str</code> with a pointer to the string
+ ** <code>"A"</code>, which mihgt be unexpected. To fix this issue,
+ ** you can use ::VL_XSTRINGIFY.
+ **
+ ** @sa ::VL_XSTRINGIFY
+ **/
+
+#define VL_STRINGIFY(x) # x
+
+/** @brief Expand and then convert the argument to a string
+ **
+ ** @param x value to be macro-expanded and converted.
+ **
+ ** This macro macro-expands the argument @a x and stringifies the
+ ** result of the expansion. For instance
+ **
+ ** @code
+ ** #define A B
+ ** char const * str = VL_STRINGIFY(A) ;
+ ** @endcode
+ **
+ ** initializes <code>str</code> with a pointer to the string
+ ** <code>"B"</code>.
+ **
+ ** @sa ::VL_STRINGIFY
+ **/
+#define VL_XSTRINGIFY(x) VL_STRINGIFY(x)
+
+/** @brief Concatenate the arguments into a lexical unit
+ **
+ ** @param x first argument to be concatenated.
+ ** @param y second argument to be concatenated.
+ **
+ ** This macro concatenates its arguments into a single lexical unit
+ ** by means of the <code>##</code> preprocessor operator. Notice that
+ ** arguments concatenated by <code>##</code> are not pre-expanded by
+ ** the C preprocessor. To macro-expand the arguments and then
+ ** concatenate them,use ::VL_XCAT.
+ **
+ ** @see ::VL_XCAT
+ **/
+#define VL_CAT(x,y) x ## y
+
+/** @brief Expand and then concatenate the arguments into a lexical unit
+ **
+ ** @param x first argument to be concatenated.
+ ** @param y second argument to be concatenated.
+ **
+ ** This macro is the same as ::VL_CAT, except that the arguments are
+ ** macro expanded before being concatenated.
+ **
+ ** @see ::VL_CAT
+ **/
+#define VL_XCAT(x,y) VL_CAT(x,y)
+
+/** @} */
+
+/** @brief Convert a boolean to "yes" or "no" strings
+ ** @param x boolean to convert.
+ ** A pointer to either the string "yes" (if @a x is true)
+ ** or the string "no".
+ ** @par Example
+ ** @code
+ ** VL_PRINTF("Is x true? %s.", VL_YESNO(x))
+ ** @endcode
+ **/
+#define VL_YESNO(x) ((x)?"yes":"no")
+
+/** ------------------------------------------------------------------
+ ** @name Heap allocation
+ ** @{ */
+
+VL_EXPORT
+void vl_set_alloc_func (void *(*malloc_func) (size_t),
+ void *(*realloc_func) (void*,size_t),
+ void *(*calloc_func) (size_t, size_t),
+ void (*free_func) (void*)) ;
+VL_INLINE void *vl_malloc (size_t n) ;
+VL_INLINE void *vl_realloc (void *ptr, size_t n) ;
+VL_INLINE void *vl_calloc (size_t n, size_t size) ;
+VL_INLINE void vl_free (void* ptr) ;
+
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Logging
+ ** @{ */
+
+/** ------------------------------------------------------------------
+ ** @brief Customizable printf function pointer type */
+typedef int(*printf_func_t) (char const *format, ...) ;
+
+/** @brief Set printf function
+ ** @param printf_func pointer to @c printf.
+ ** Let @c print_func be NULL to disable printf.
+ **/
+VL_EXPORT void vl_set_printf_func (printf_func_t printf_func) ;
+
+/** @def VL_PRINTF
+ ** @brief Call user-customizable @c printf function
+ **
+ ** The function calls the user customizable @c printf.
+ **/
+#define VL_PRINTF (*vl_printf_func)
+
+/** @def VL_PRINT
+ ** @brief Same as ::VL_PRINTF (legacy code)
+ **/
+#define VL_PRINT (*vl_printf_func)
+
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Error handling
+ ** @{ */
+
+/** @brief The number of the last error */
+extern VL_EXPORT int vl_err_no ;
+
+/** @brief The maximum length of an error description. */
+#define VL_ERR_MSG_LEN 1024
+
+/** @brief The description of the last error. */
+extern VL_EXPORT char vl_err_msg [VL_ERR_MSG_LEN + 1] ;
+
+#define VL_ERR_OK 0 /**< No error */
+#define VL_ERR_OVERFLOW 1 /**< Buffer overflow error */
+#define VL_ERR_ALLOC 2 /**< Resource allocation error */
+#define VL_ERR_BAD_ARG 3 /**< Bad argument or illegal data error */
+#define VL_ERR_IO 4 /**< Input/output error */
+#define VL_ERR_EOF 5 /**< End-of-file or end-of-sequence error */
+#define VL_ERR_NO_MORE 5 /**< End-of-sequence @deprecated */
+
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Common operations
+ ** @{ */
+
+/** @brief Min operation
+ ** @param x value
+ ** @param y value
+ ** @return the minimum of @a x and @a y.
+ **/
+#define VL_MIN(x,y) (((x)<(y))?(x):(y))
+
+/** @brief Max operation
+ ** @param x value.
+ ** @param y value.
+ ** @return the maximum of @a x and @a y.
+ **/
+#define VL_MAX(x,y) (((x)>(y))?(x):(y))
+
+/** @brief Signed left shift operation
+ **
+ ** The macro is equivalent to the builtin @c << operator, but it
+ ** supports negative shifts too.
+ **
+ ** @param x value.
+ ** @param n number of shift positions.
+ ** @return @c x << n .
+ **/
+#define VL_SHIFT_LEFT(x,n) (((n)>=0)?((x)<<(n)):((x)>>-(n)))
+/* @} */
+
+/** --------------------------------------------------------------- */
+VL_EXPORT
+char const * vl_get_version_string () ;
+
+VL_EXPORT
+void vl_print_info () ;
+
+/** ------------------------------------------------------------------
+ ** @name Measuring time
+ ** @{
+ **/
+VL_EXPORT void vl_tic() ;
+VL_EXPORT double vl_toc() ;
+/** @} */
+
+extern VL_EXPORT int (*vl_printf_func) (char const * format, ...) ;
+extern VL_EXPORT void *(*vl_malloc_func) (size_t) ;
+extern VL_EXPORT void *(*vl_realloc_func) (void*,size_t) ;
+extern VL_EXPORT void *(*vl_calloc_func) (size_t, size_t) ;
+extern VL_EXPORT void (*vl_free_func) (void*) ;
+
+VL_INLINE
+void* vl_malloc (size_t n)
+{
+ return (*vl_malloc_func)(n) ;
+}
+
+VL_INLINE
+void* vl_realloc (void* ptr, size_t n)
+{
+ return (*vl_realloc_func)(ptr, n) ;
+}
+
+VL_INLINE
+void* vl_calloc (size_t n, size_t size)
+{
+ return (*vl_calloc_func)(n, size) ;
+}
+
+VL_INLINE
+void vl_free (void *ptr)
+{
+ (*vl_free_func)(ptr) ;
+}
+
+/* VL_GENERIC_H */
+#endif
diff --git a/vl/host.c b/vl/host.c
new file mode 100644
index 0000000..18f4ec2
--- /dev/null
+++ b/vl/host.c
@@ -0,0 +1,509 @@
+/** @file host.c
+ ** @author Andrea Vedaldi
+ ** @brief Host - Definition
+ **/
+
+/* AUTORIGHTS
+ Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+ This file is part of VLFeat, available in the terms of the GNU
+ General Public License version 2.
+ */
+
+#include "host.h"
+#include "generic.h"
+
+/**
+ @file host.h
+
+ This module provides functionalities to identify the
+ host operating system, C compiler,
+ and CPU architecture. It also provides a few features
+ to abstract from such details.
+
+ - @ref host-os "Host operating system"
+ - @ref host-compiler "Host compiler"
+ - @ref host-compiler-data-model "Data models"
+ - @ref host-compiler-other "Oter compiler specific features"
+ - @ref host-arch "Host CPU architecture"
+ - @ref host-arch-endianness "Endianness"
+
+ @see http://predef.sourceforge.net/index.php
+
+ @section host-os Host operating system
+
+ The module defines a symbol to identify the host operating system:
+ ::VL_OS_WIN for Windows, ::VL_OS_LINUX for Linux, ::VL_OS_MACOSX for
+ Mac OS X, and so on.
+
+ @section host-compiler Host compiler
+
+ The module defines a symbol to identify the host compiler:
+ ::VL_COMPILER_MSC for Microsoft Visual C++, ::VL_COMPILER_GNUC for GNU C,
+ and so on. The (integer) value of such symbols
+ corresponds the version of the compiler.
+
+ The module defines a symbol to identify the data model of the compiler:
+ ::VL_COMPILER_ILP32, ::VL_COMPILER_LP32, or ::VL_COMPILER_LP64 (see
+ Sect. @ref host-compiler-data-model).
+ For convenience, it also defines a number of atomic types of prescribed
+ width (::vl_int8, ::vl_int16, ::vl_int32, etc.).
+
+ @remark While some of such functionalities are provided by the
+ standard header @c stdint.h, the latter is not supported
+ by all platforms.
+
+ @subsection host-compiler-data-model Data models
+
+ The C language defines a number of atomic data types (such as @c char, @c short,
+ @c int and so on). The number of bits (width)
+ used to represent each data type depends on the compiler data model.
+ The following table summarizes the relevant conventions:
+
+ <table><caption><b>Compiler data models.</b> The table shows
+ how many bits are allocated to each atomic data type according to
+ each model.</caption>
+ <tr style="font-weight:bold;">
+ <td>Data model</td>
+ <td><code>short</code></td>
+ <td><code>int</code></td>
+ <td><code>long</code></td>
+ <td><code>long long</code></td>
+ <td><code>void*</code></td>
+ <td>Compiler</td>
+ </tr>
+ <tr>
+ <td>ILP32</td>
+ <td style="background-color:#ffa;">16</td>
+ <td style="background-color:#afa;">32</td>
+ <td style="background-color:#afa;">32</td>
+ <td >64</td>
+ <td style="background-color:#afa;">32</td>
+ <td>common 32 bit architectures</td>
+ </tr>
+ <tr>
+ <td>LP64</td>
+ <td style="background-color:#ffa;">16</td>
+ <td style="background-color:#afa;">32</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>UNIX-64 (Linux, Mac OS X)</td>
+ </tr>
+ <tr>
+ <td>ILP64</td>
+ <td style="background-color:#ffa;">16</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>Alpha, Cray</td>
+ </tr>
+ <tr>
+ <td>SLIP64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td>64</td>
+ <td></td>
+ </tr>
+ <tr>
+ <td>LLP64</td>
+ <td style="background-color:#ffa;">16</td>
+ <td style="background-color:#afa;">32</td>
+ <td style="background-color:#afa;">32</td>
+ <td>64</td>
+ <td>64</td>
+ <td>Windows-64</td>
+ </tr>
+ </table>
+
+ @subsection host-compiler-other Other compiler-specific features
+
+ The module provides the macro ::VL_EXPORT to declare symbols exported
+ from
+ the library and the macro ::VL_INLINE to declare inline functions.
+ Such features are not part of the C89 standard, and change
+ depending on the compiler.
+
+ @par "Example:"
+ The following header file declares a function @c f that
+ should be visible from outside the library.
+ @code
+ #include <vl/generic.h>
+ VL_EXPORT void f () ;
+ VL_EXPORT int i ;
+ @endcode
+ Notice that the macro ::VL_EXPORT needs not to be included again
+ when the function is defined.
+
+ @par "Example:"
+ The following header file declares an inline function @c f:
+ @code
+ #include <vl/generic.h>
+ VL_INLINE int f() ;
+
+ VL_INLINE int f() { return 1 ; }
+ @endcode
+ Here the first instruction defines the function @c f, where the second
+ declares it. Notice that since this is an inline function, its definition
+ must be found in the header file rather than in an implementation file.
+ Notice also that definition and declaration can be merged.
+
+ These macros translate according to the following tables:
+
+ <table style="font-size:70%;">
+ <caption>Macros for exporting library symbols</caption>
+ <tr>
+ <td>Platform</td>
+ <td>Macro name</td>
+ <td>Value when building the library</td>
+ <td>Value when importing the library</td>
+ </tr>
+ <tr>
+ <td>Unix/GCC</td>
+ <td>::VL_EXPORT</td>
+ <td>empty (assumes <c>-visibility=hidden</c> GCC option)</td>
+ <td><c>__attribute__((visibility ("default")))</c></td>
+ </tr>
+ <tr>
+ <td>Win/Visual C++</td>
+ <td>::VL_EXPORT</td>
+ <td>@c __declspec(dllexport)</td>
+ <td>@c __declspec(dllimport)</td>
+ </tr>
+ </table>
+
+ <table style="font-size:70%;">
+ <caption>Macros for declaring inline functions</caption>
+ <tr>
+ <td>Platform</td>
+ <td>Macro name</td>
+ <td>Value</td>
+ </tr>
+ <tr>
+ <td>Unix/GCC</td>
+ <td>::VL_INLINE</td>
+ <td>static inline</td>
+ </tr>
+ <tr>
+ <td>Win/Visual C++</td>
+ <td>::VL_INLINE</td>
+ <td>static __inline</td>
+ </tr>
+ </table>
+
+ @section host-arch Host CPU architecture
+
+ The module defines a symbol to identify the host CPU architecture:
+ ::VL_ARCH_IX86 for Intel x86, ::VL_ARCH_IA64 for Intel 64, and so on.
+
+ @subsection host-arch-endianness Endianness
+
+ The module defines a symbol to identify the host CPU endianness:
+ ::VL_ARCH_BIG_ENDIAN for big endian and ::VL_ARCH_LITTLE_ENDIAN for
+ little endian. The functions
+ ::vl_swap_host_big_endianness_8(), ::vl_swap_host_big_endianness_4(),
+ ::vl_swap_host_big_endianness_2() to change the endianness
+ of data (from/to host and network order).
+
+ Recall that <em>endianness</em> concerns the way multi-byte data types
+ (such as 16, 32 and 64 bits
+ integers) are stored into the addressable memory.
+ All CPUs uses a contiguous address range to store atomic data types
+ (e.g. a 16-bit integer could
+ be assigned to the addresses <c>0x10001</c> and <c>0x10002</c>), but
+ the order may differ.
+
+ - The convention is <em>big endian</em>, or in <em>network order</em>,
+ if the most significant byte of the multi-byte data types is assigned to
+ the smaller memory address. This is the convention used for instance
+ by the PPC architecture.
+
+ - The convention is <em>little endian</em> if the least significant
+ byte is assigned to the smaller memory address. This is the convention
+ used for instance by the x86 architecture.
+
+ @remark The names &ldquo;big endian&rdquo; and &ldquo;little endian&rdquo; are
+ a little confusing. &ldquo;Big endian&rdquo; means &ldquo;big endian first&rdquo;, i.e.
+ the address of the most significant byte comes first. Similarly,
+ &ldquo;little endian&rdquo; means &ldquo;little endian first&rdquo;,
+ in the sense that the address of the least significant byte comes first.
+
+ Endianness is a concern when data is either exchanged with processors
+ that use different conventions,
+ transmitted over a network, or stored
+ to a file. For the latter two cases, one usually saves data in
+ big endian (network) order regardless of the host CPU.
+
+ **/
+
+/** @def VL_OS_LINUX
+ ** @brief Defined if the host operating system is Linux.
+ **/
+
+/** @def VL_OS_MACOSX
+ ** @brief Defined if the host operating system is Mac OS X.
+ **/
+
+/** @def VL_OS_WIN
+ ** @brief Defined if the host operating system is Windows (32 or 64)
+ **/
+
+/** @def VL_OS_WIN64
+ ** @brief Defined if the host operating system is Windows-64.
+ **/
+
+/** @def VL_COMPILER_GNUC
+ ** @brief Defined if the host compiler is GNU C.
+ **
+ ** This macro is defined if the compiler is GNUC.
+ ** Its value is calculated as
+ ** @code
+ ** 10000 * MAJOR + 100 * MINOR + PATCHLEVEL
+ ** @endcode
+ ** @see @ref host-compiler
+ **/
+
+/** @def VL_COMPILER_MSC
+ ** @brief Defined if the host compiler is Microsoft Visual C++.
+ ** @see @ref host-compiler
+ **/
+
+/** @def VL_COMPILER_LLP64
+ ** @brief Defined if the host compiler data model is LLP64.
+ ** @see @ref host-compiler-data-model
+ **/
+
+/** @def VL_COMPILER_LP64
+ ** @brief Defined if the host compiler data model is LP64.
+ ** @see @ref host-compiler-data-model
+ **/
+
+/** @def VL_COMPILER_ILP32
+ ** @brief Defined if the host compiler data model is ILP32.
+ ** @see @ref host-compiler-data-model
+ **/
+
+/** @def VL_ARCH_IX86
+ ** @brief Defined if the host CPU is of the Intel x86 family.
+ ** @see @ref host-arch
+ **/
+
+/** @def VL_ARCH_IA64
+ ** @brief Defined if the host CPU is of the Intel Architecture-64 family.
+ ** @see @ref host-arch
+ **/
+
+/** @def VL_ARCH_LITTLE_ENDIAN
+ ** @brief Defined if the host CPU is little endian
+ ** @see @ref host-arch-endianness
+ **/
+
+/** @def VL_ARCH_BIG_ENDIAN
+ ** @brief Defined if the host CPU is big endian
+ ** @see @ref host-arch-endianness
+ **/
+
+/** @def VL_INLINE
+ ** @brief Adds appropriate inline function qualifier
+ ** @see @ref host-compiler-others
+ **/
+
+/** @def VL_EXPORT
+ ** @brief Declares a DLL exported symbol
+ ** @see @ref host-compiler-others
+ **/
+
+/** --------------------------------------------------------------- */
+
+#if defined(VL_ARCH_IX86) || defined(VL_ARCH_IA64)
+#define HAS_CPUID
+#endif
+
+#if defined(HAS_CPUID) & defined(VL_COMPILER_MSC)
+#include <intrin.h>
+VL_INLINE
+void _vl_cpuid (vl_int32* info, int function)
+{
+ __cpuid(info, function) ;
+}
+#endif
+
+#if defined(HAS_CPUID) & defined(VL_COMPILER_GNUC)
+VL_INLINE
+void _vl_cpuid (vl_int32* info, int function)
+{
+ /* this version is compatible with -fPIC */
+ __asm__ __volatile__
+ ("pushl %%ebx \n\t" /* save %ebx */
+ "cpuid \n\t"
+ "movl %%ebx, %1 \n\t" /* save what cpuid just put in %ebx */
+ "popl %%ebx \n\t" /* restore the old %ebx */
+ : "=a"(info[0]), "=r"(info[1]), "=c"(info[2]), "=d"(info[3])
+ : "a"(function)
+ : "cc") ;
+}
+#endif
+
+#ifdef HAS_CPUID
+struct x86cpu_
+{
+ char vendor_string [0x20] ;
+ vl_bool has_sse42 ;
+ vl_bool has_sse41 ;
+ vl_bool has_sse3 ;
+ vl_bool has_sse2 ;
+ vl_bool has_sse ;
+ vl_bool has_mmx ;
+} x86cpu ;
+
+vl_bool x86cpu_initialized = 0 ;
+
+void _vl_x86cpu_init ()
+{
+ vl_int32 info [4] ;
+ int max_func = 0 ;
+ _vl_cpuid(info, 0) ;
+ max_func = info[0] ;
+ *((vl_int32*)x86cpu.vendor_string+0) = info[1] ;
+ *((vl_int32*)x86cpu.vendor_string+1) = info[3] ;
+ *((vl_int32*)x86cpu.vendor_string+2) = info[2] ;
+
+ if (max_func >= 1) {
+ _vl_cpuid(info, 1) ;
+ x86cpu.has_mmx = info[3] & (1 << 23) ;
+ x86cpu.has_sse = info[3] & (1 << 25) ;
+ x86cpu.has_sse2 = info[3] & (1 << 26) ;
+ x86cpu.has_sse3 = info[2] & (1 << 0) ;
+ x86cpu.has_sse41 = info[2] & (1 << 19) ;
+ x86cpu.has_sse42 = info[2] & (1 << 20) ;
+ }
+}
+
+VL_INLINE
+struct x86cpu_ const* _vl_x86cpu_get()
+{
+ if (!x86cpu_initialized) _vl_x86cpu_init() ;
+ return &x86cpu ;
+}
+#endif
+
+/** --------------------------------------------------------------- */
+
+vl_bool simd_enabled = 1 ;
+
+/** @brief Enalbe/disable usage of SIMD instructions
+ ** @param x set to @c true to enable SIMD instructions.
+ **
+ ** Notice that usage of SIMD instructions may be prevented due
+ ** to lack of CPU support and data alignment issues.
+ **
+ ** @see ::vl_cpu_has_sse2(), ::vl_cpu_has_sse3(), etc.
+ **/
+void vl_set_simd_enabled (vl_bool x)
+{
+ simd_enabled = x ;
+}
+
+/** @brief Are SIMD instructons enabled?
+ ** @return @c true is SIMD is enabled.
+ **/
+vl_bool vl_get_simd_enabled ()
+{
+ return simd_enabled ;
+}
+
+/** @brief Check for SSE3 instruction set
+ ** @return @c true if SSE3 is present.
+ **/
+vl_bool vl_cpu_has_sse3 ()
+{
+#ifdef HAS_CPUID
+ return _vl_x86cpu_get()->has_sse3 ;
+#else
+ return 0 ;
+#endif
+}
+
+/** @brief Check for SSE2 instruction set
+ ** @return @c true if SSE2 is present.
+ **/
+vl_bool vl_cpu_has_sse2 ()
+{
+#ifdef HAS_CPUID
+ return _vl_x86cpu_get()->has_sse2 ;
+#else
+ return 0 ;
+#endif
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Print host information
+ **/
+
+void vl_print_host_info ()
+{
+ char const *arch = 0, *endian = 0, *comp = 0, *dm = 0 ;
+ int compver ;
+#ifdef VL_ARCH_IA6
+ arch = "IA64" ;
+#endif
+#ifdef VL_ARCH_IX86
+ arch = "IX86" ;
+#endif
+#ifdef VL_ARCH_PPC
+ arch = "PPC" ;
+#endif
+
+#ifdef VL_ARCH_BIN_ENDIAN
+ endian = "big endian" ;
+#endif
+#ifdef VL_ARCH_LITTLE_ENDIAN
+ endian = "little endian" ;
+#endif
+
+#ifdef VL_COMPILER_MSC
+ comp = "Microsoft Visual C++" ;
+ compver = VL_COMPILER_MSC ;
+#endif
+#ifdef VL_COMPILER_GNUC
+ comp = "GNU C" ;
+ compver = VL_COMPILER_GNUC ;
+#endif
+
+#ifdef VL_COMPILER_LP64
+ dm = "LP64" ;
+#endif
+#ifdef VL_COMPILER_LLP64
+ dm = "LP64" ;
+#endif
+#ifdef VL_COMPILER_ILP32
+ dm = "ILP32" ;
+#endif
+
+#define YESNO(x) ((x)?"yes":"no")
+
+ VL_PRINTF("Host: Compiler: %s %d\n", comp, compver) ;
+ VL_PRINTF(" Compiler data model: %s\n", dm) ;
+ VL_PRINTF(" CPU architecture: %s\n", arch) ;
+ VL_PRINTF(" CPU endianness: %s\n", endian) ;
+
+#ifdef HAS_CPUID
+ {
+ struct x86cpu_ const* c = _vl_x86cpu_get() ;
+ VL_PRINTF(" CPU vendor string: %s\n", c->vendor_string) ;
+ VL_PRINTF(" CPU has MMX: %s\n", YESNO(c->has_mmx)) ;
+ VL_PRINTF(" CPU has SSE: %s\n", YESNO(c->has_sse)) ;
+ VL_PRINTF(" CPU has SSE2: %s\n", YESNO(c->has_sse2)) ;
+ VL_PRINTF(" CPU has SSE3: %s\n", YESNO(c->has_sse3)) ;
+ VL_PRINTF(" CPU has SSE4.1: %s\n", YESNO(c->has_sse41)) ;
+ VL_PRINTF(" CPU has SSE4.2: %s\n", YESNO(c->has_sse42)) ;
+ VL_PRINTF("VLFeat uses SIMD: %s\n", YESNO(vl_get_simd_enabled())) ;
+ }
+#endif
+}
+
+
+
diff --git a/vl/host.h b/vl/host.h
new file mode 100644
index 0000000..92580e9
--- /dev/null
+++ b/vl/host.h
@@ -0,0 +1,367 @@
+/** @file host.h
+ ** @author Andrea Vedaldi
+ ** @brief Host
+ **/
+
+#ifndef VL_HOST_H
+#define VL_HOST_H
+
+/*
+ The following macros identify the host OS, architecture and compiler.
+ They are derived from http://predef.sourceforge.net/
+ */
+
+/** @name Identifying the host operating system
+ ** @{ */
+#if defined(linux) || \
+ defined(__linux) || \
+ defined(__linux__) || \
+ defined(__DOXYGEN__)
+#define VL_OS_LINUX 1
+#endif
+
+#if (defined(__APPLE__) & defined(__MACH__)) || \
+ defined(__DOXYGEN__)
+#define VL_OS_MACOSX 1
+#endif
+
+#if defined(__WIN32__) || \
+ defined(_WIN32) || \
+ defined(__DOXYGEN__)
+#define VL_OS_WIN 1
+#endif
+
+#if defined(_WIN64) || \
+ defined(__DOXYGEN__)
+#define VL_OS_WIN64 1
+#endif
+/** @} */
+
+/** @name Identifying the host compiler
+ ** @{ */
+#if defined(__GNUC__) || defined(__DOXYGEN__)
+# if defined(__GNUC_PATCHLEVEL__)
+# define VL_COMPILER_GNUC (__GNUC__ * 10000 \
++ __GNUC_MINOR__ * 100 \
++ __GNUC_PATCHLEVEL__)
+# else
+# define VL_COMPILER_GNUC (__GNUC__ * 10000 \
++ __GNUC_MINOR__ * 100)
+# endif
+#endif
+
+#if defined(_MSC_VER) || defined(__DOXYGEN__)
+#define VL_COMPILER_MSC _MSC_VER
+#endif
+
+#if defined(__LCC__) || defined(__DOXYGEN__)
+#warning "LCC support is experimental!"
+#define VL_COMPILER_LCC 1
+#endif
+
+/** @} */
+
+/** @name Identifying the host CPU architecture
+ ** @{ */
+#if defined(i386) || \
+ defined(__i386__) || \
+ defined(__DOXYGEN__)
+#define VL_ARCH_IX86 300
+#elif defined(__i486__)
+#define VL_ARCH_IX86 400
+#elif defined(__i586__)
+#define VL_ARCH_IX86 500
+#elif defined(__i686__)
+#define VL_ARCH_IX86 600
+#elif defined(_M_IX86)
+#define VL_ARCH_IX86 _M_IX86
+#endif
+
+#if defined(__ia64__) || \
+ defined(_IA64) || \
+ defined(__IA64) || \
+ defined(__ia64) || \
+ defined(_M_IA64) || \
+ defined(__DOXYGEN__)
+#define VL_ARCH_IA64
+#endif
+/** @} */
+
+/** @name Identifying the host data model
+ ** @{ */
+#if defined(__LLP64__) || \
+ defined(__LLP64) || \
+ defined(__LLP64) || \
+ (defined(VL_COMPILER_MSC) & defined(VL_OS_WIN64)) || \
+ (defined(VL_COMPILER_LCC) & defined(VL_OS_WIN64)) || \
+ defined(__DOXYGEN__)
+#define VL_COMPILER_LLP64
+#endif
+
+#if defined(__LP64__) || \
+ defined(__LP64) || \
+ defined(__LP64) || \
+ (defined(VL_OS_MACOSX) & defined(VL_ARCH_IA64)) || \
+ defined(__DOXYGEN__)
+#define VL_COMPILER_LP64
+#endif
+
+#if (!defined(VL_COMPILER_LLP64) & !defined(VL_COMPILER_LP64)) || \
+ defined(__DOXYGEN__)
+#define VL_COMPILER_ILP32
+#endif
+/** @} */
+
+/** @name Identifying the host endianness
+ ** @{ */
+#if defined(__LITTLE_ENDIAN__) || \
+ defined(VL_ARCH_IX86) || \
+ defined(VL_ARCH_IA64) || \
+ defined(__DOXYGEN__)
+#define VL_ARCH_LITTLE_ENDIAN
+#endif
+
+#if defined(__DOXYGEN__) || \
+ !defined(VL_ARCH_LITTLE_ENDIAN)
+#define VL_ARCH_BIG_ENDIAN
+#endif
+/** @} */
+
+#if defined(VL_COMPILER_MSC)
+#define VL_INLINE static __inline
+#define snprintf _snprintf
+#define isnan _isnan
+#ifdef VL_BUILD_DLL
+#define VL_EXPORT __declspec(dllexport)
+#else
+#define VL_EXPORT //__declspec(dllimport)
+#endif
+#endif
+
+#if defined(VL_COMPILER_LCC)
+#define VL_INLINE static __inline
+#define snprintf _snprintf
+#define isnan _isnan
+#ifdef VL_BUILD_DLL
+#define VL_EXPORT __declspec(dllexport)
+#else
+#define VL_EXPORT
+#endif
+#endif
+
+#if defined(VL_COMPILER_GNUC) || \
+ defined(__DOXYGEN__)
+#define VL_INLINE static __inline__
+#ifdef VL_BUILD_DLL
+#define VL_EXPORT __attribute__((visibility ("default")))
+#else
+#define VL_EXPORT
+#endif
+#endif
+
+/** ------------------------------------------------------------------
+ ** @name Atomic data types
+ ** @{
+ **/
+
+#define VL_TRUE 1 /**< @brief @c true (1) constant */
+#define VL_FALSE 0 /**< @brief @c false (0) constant */
+
+#if defined(VL_COMPILER_LP64) || defined(VL_COMPILER_LLP64)
+typedef long long vl_int64 ; /**< @brief Signed 64-bit integer. */
+typedef int vl_int32 ; /**< @brief Signed 32-bit integer. */
+typedef short vl_int16 ; /**< @brief Signed 16-bit integer. */
+typedef char vl_int8 ; /**< @brief Signed 8-bit integer. */
+
+typedef long long unsigned vl_uint64 ; /**< @brief Unsigned 64-bit integer. */
+typedef int unsigned vl_uint32 ; /**< @brief Unsigned 32-bit integer. */
+typedef short unsigned vl_uint16 ; /**< @brief Unsigned 16-bit integer. */
+typedef char unsigned vl_uint8 ; /**< @brief Unsigned 8-bit integer. */
+
+typedef int vl_int ; /**< @brief Same as @c int. */
+typedef unsigned int vl_uint ; /**< @brief Same as <code>unsigned int</code>. */
+
+typedef int vl_bool ; /**< @brief Boolean. */
+typedef vl_int64 vl_intptr ; /**< @brief Integer holding a pointer. */
+typedef vl_uint64 vl_uintptr ; /**< @brief Unsigned integer holding a pointer. */
+#endif
+
+#if defined(VL_COMPILER_ILP32)
+
+#ifdef VL_COMPILER_MSC
+typedef __int64 vl_int64 ;
+#else
+typedef long long vl_int64 ;
+#endif
+
+typedef int vl_int32 ;
+typedef short vl_int16 ;
+typedef char vl_int8 ;
+
+#ifdef VL_COMPILER_MSC
+typedef __int64 unsigned vl_uint64 ;
+#else
+typedef long long unsigned vl_uint64 ;
+#endif
+typedef int unsigned vl_uint32 ;
+typedef short unsigned vl_uint16 ;
+typedef char unsigned vl_uint8 ;
+
+typedef int vl_int ;
+typedef unsigned int vl_uint ;
+
+typedef int vl_bool ;
+typedef vl_int32 vl_intptr ;
+typedef vl_uint32 vl_uintptr ;
+#endif
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Printing the atomic data types
+ ** @{ */
+
+
+/** @def VL_FL_INT64
+ ** @brief @c prinf length flag for ::vl_int64 and ::vl_uint64.
+ **/
+
+/** @def VL_FL_INT32
+ ** @brief @c prinf length flag for ::vl_int32 and ::vl_uint32.
+ **/
+
+/** @def VL_FL_INT16
+ ** @brief @c prinf length flag for ::vl_int16 and ::vl_uint16.
+ **/
+
+/** @def VL_FL_INT8
+ ** @brief @c prinf length flag for ::vl_int8 and ::vl_uint8.
+ **/
+#ifdef VL_COMPILER_MSC
+#define VL_FL_INT64 "I64"
+#else
+#define VL_FL_INT64 "ll"
+#endif
+#define VL_FL_INT32 ""
+#define VL_FL_INT16 "h"
+#define VL_FL_INT8 "hh"
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Atomic data types limits
+ ** @{ */
+
+/** @brief Largest integer (math constant) */
+#define VL_BIG_INT 0x7FFFFFFFL
+
+/** @brief Smallest integer (math constant) */
+#define VL_SMALL_INT (- VL_BIG_INT - 1)
+
+/** @brief Largest unsigned integer (math constant) */
+#define VL_BIG_UINT 0xFFFFFFFFUL
+
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Endianness detection and conversion
+ ** @{
+ **/
+VL_INLINE void vl_swap_host_big_endianness_8 (void *dst, void* src) ;
+VL_INLINE void vl_swap_host_big_endianness_4 (void *dst, void* src) ;
+VL_INLINE void vl_swap_host_big_endianness_2 (void *dst, void* src) ;
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @name Obtaining host info at run time
+ ** @{ */
+VL_EXPORT void vl_print_host_info() ;
+VL_EXPORT vl_bool vl_cpu_has_sse3 () ;
+VL_EXPORT vl_bool vl_cpu_has_sse2 () ;
+/** @} */
+
+VL_EXPORT void vl_set_simd_enabled (vl_bool x) ;
+VL_EXPORT vl_bool vl_get_simd_enabled() ;
+
+/** ------------------------------------------------------------------
+ ** @brief Host <-> big endian transformation for 8-bytes value
+ **
+ ** @param dst destination 8-byte buffer.
+ ** @param src source 8-byte bufffer.
+ ** @see @ref host-arch-endianness.
+ **/
+
+VL_INLINE void
+vl_swap_host_big_endianness_8 (void *dst, void* src)
+{
+ char *dst_ = (char*) dst ;
+ char *src_ = (char*) src ;
+#if defined(VL_ARCH_BIG_ENDIAN)
+ dst_ [0] = src_ [0] ;
+ dst_ [1] = src_ [1] ;
+ dst_ [2] = src_ [2] ;
+ dst_ [3] = src_ [3] ;
+ dst_ [4] = src_ [4] ;
+ dst_ [5] = src_ [5] ;
+ dst_ [6] = src_ [6] ;
+ dst_ [7] = src_ [7] ;
+#else
+ dst_ [0] = src_ [7] ;
+ dst_ [1] = src_ [6] ;
+ dst_ [2] = src_ [5] ;
+ dst_ [3] = src_ [4] ;
+ dst_ [4] = src_ [3] ;
+ dst_ [5] = src_ [2] ;
+ dst_ [6] = src_ [1] ;
+ dst_ [7] = src_ [0] ;
+#endif
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Host <-> big endian transformation for 4-bytes value
+ **
+ ** @param dst destination 4-byte buffer.
+ ** @param src source 4-byte bufffer.
+ ** @sa @ref host-arch-endianness.
+ **/
+
+VL_INLINE void
+vl_swap_host_big_endianness_4 (void *dst, void* src)
+{
+ char *dst_ = (char*) dst ;
+ char *src_ = (char*) src ;
+#if defined(VL_ARCH_BIG_ENDIAN)
+ dst_ [0] = src_ [0] ;
+ dst_ [1] = src_ [1] ;
+ dst_ [2] = src_ [2] ;
+ dst_ [3] = src_ [3] ;
+#else
+ dst_ [0] = src_ [3] ;
+ dst_ [1] = src_ [2] ;
+ dst_ [2] = src_ [1] ;
+ dst_ [3] = src_ [0] ;
+#endif
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Host <-> big endian transformation for 2-bytes value
+ **
+ ** @param dst destination 2-byte buffer.
+ ** @param src source 2-byte bufffer.
+ ** @see @ref host-arch-endianness.
+ **/
+
+VL_INLINE void
+vl_swap_host_big_endianness_2 (void *dst, void* src)
+{
+ char *dst_ = (char*) dst ;
+ char *src_ = (char*) src ;
+#if defined(VL_ARCH_BIG_ENDIAN)
+ dst_ [0] = src_ [0] ;
+ dst_ [1] = src_ [1] ;
+#else
+ dst_ [0] = src_ [1] ;
+ dst_ [1] = src_ [0] ;
+#endif
+}
+
+/* VL_HOST_H */
+#endif
diff --git a/vl/imop.c b/vl/imop.c
new file mode 100644
index 0000000..69899ef
--- /dev/null
+++ b/vl/imop.c
@@ -0,0 +1,85 @@
+#if WIN32
+#pragma warning(disable:4244)
+#endif
+
+/** @file imop.c
+ ** @author Andrea Vedaldi
+ ** @brief Image operations - Definition
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#include "imop.h"
+
+#include <assert.h>
+#include <string.h>
+
+/** @fn ::vl_convtransp_f(
+ ** float*,float const*,float const*,int,int,int,int)
+ **
+ ** @brief Convolve along columns and take transpose
+ **
+ ** The function convolves the columns of the matrix @a src by the
+ ** kernel @a filt and writes the transpose of the
+ ** result to the buffer @a dst.
+ **
+ ** @remark Here `columns' correspond to the fastest varying
+ ** index. Often image are stored in row-major order, so the `columns'
+ ** are actually the `rows' of the image.
+ **
+ ** @remark To convolve an image by a 2-D separable filter just call
+ ** this function twice.
+ **
+ ** @remark The function can be used to process N-dimensional arrays
+ ** by stacking them properly.
+ **
+ ** @param dst output image buffer.
+ ** @param src input image buffer.
+ ** @param filt filter buffer.
+ ** @param width width of the image.
+ ** @param height height of the image.
+ ** @param filt_width filter half width (total width = 2 @c filt_width + 1).
+ ** @param mode mode flags.
+ **/
+
+/** @fn ::vl_imsmooth_f(
+ ** float*,float*,float const*,int,int,double)
+ **
+ ** @brief Smooth image by Gaussian kernel
+ **
+ ** The function convolves the image @a src by a Gaussian kernel of
+ ** variance @a sigma and write the result to the buffer @a dst. The
+ ** functions also need a scratch buffer @a dst of the same size of
+ ** the buffers @a src and @a dst.
+ **
+ ** The various functions differ by the data type of the image pixels
+ ** that they process.
+ **
+ ** @param dst output image buffer.
+ ** @param temp scratch image buffer.
+ ** @param src input image buffer.
+ ** @param width width of the buffers.
+ ** @param height height of the buffers.
+ ** @param sigma standard deviation of the Gaussian kernel.
+ **/
+
+#define PIX float /**< pixel type. @internal */
+#define FLT float /**< float type. @internal */
+#define SFX f /**< suffix. @internal */
+#include "imop.tc"
+#undef PIX
+#undef SFX
+#undef FLT
+
+#define PIX double
+#define FLT double
+#define SFX d
+#include "imop.tc"
+#undef PIX
+#undef SFX
+#undef FLT
diff --git a/vl/imop.h b/vl/imop.h
new file mode 100644
index 0000000..6d59414
--- /dev/null
+++ b/vl/imop.h
@@ -0,0 +1,62 @@
+/** @file imop.h
+ ** @author Andrea Vedaldi
+ ** @brief Image operations
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#ifndef VL_IMOP
+#define VL_IMOP
+
+#include "generic.h"
+
+/** @name Convolution padding modes
+ ** @{
+ **/
+#define VL_CONV_ZERO 0 /**< Zero padding */
+#define VL_CONV_CIRC 1 /**< Circular convolution */
+#define VL_CONV_CONT 2 /**< Pad by continuity */
+/* @} */
+
+/** @name Convolution
+ ** @{
+ **/
+VL_EXPORT
+void vl_convtransp_f(float *dst,
+ float const *src,
+ float const *filt,
+ int width, int height, int filt_width,
+ int mode) ;
+
+VL_EXPORT
+void vl_convtransp_d(double *dst,
+ double const *src,
+ double const *filt,
+ int width, int height, int filt_width,
+ int mode) ;
+/* @} */
+
+
+/** @name Image Smoothing
+ ** @{
+ **/
+VL_EXPORT
+void vl_imsmooth_f(float *dst,
+ float *temp,
+ float const *src,
+ int width, int height, double sigma) ;
+
+VL_EXPORT
+void vl_imsmooth_d(double *dst,
+ double *temp,
+ double const *src,
+ int width, int height, double sigma) ;
+/*@}*/
+
+/* VL_IMOP */
+#endif
diff --git a/vl/imop.tc b/vl/imop.tc
new file mode 100644
index 0000000..aa9d496
--- /dev/null
+++ b/vl/imop.tc
@@ -0,0 +1,155 @@
+/** @file imop.tc
+ ** @author Andrea Vedaldi
+ ** @brief Image operations - Definition Template
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#include "mathop.h"
+
+#include <math.h>
+#include <stdlib.h>
+
+#undef VL_IMSMOOTH
+#undef VL_IMSMOOTH_
+#undef VL_CONVTRANSP
+#undef VL_CONVTRANSP_
+
+#define VL_IMSMOOTH(SFX) vl_imsmooth_ ## SFX
+#define VL_IMSMOOTH_(SFX) VL_IMSMOOTH(SFX)
+#define VL_CONVTRANSP(SFX) vl_convtransp_ ## SFX
+#define VL_CONVTRANSP_(SFX) VL_CONVTRANSP(SFX)
+
+VL_EXPORT
+void
+VL_CONVTRANSP_(SFX) (PIX *dst,
+ PIX const *src,
+ PIX const *filt,
+ int width, int height, int filt_width,
+ int mode)
+{
+ int i, j ;
+
+ /* Convolve along the first dimension. Also, circularly swap */
+ /* dimension in the output, making the first dimension the last */
+ /* */
+ /* src is width x height */
+ /* dst is height x width */
+ /* filt is 2 * filt_width + 1 */
+ /* */
+
+ switch (mode) {
+ case VL_CONV_CONT :
+ for(j = 0 ; j < height ; ++j) {
+ for(i = 0 ; i < width ; ++i) {
+ PIX acc = 0.0 ;
+ PIX const *g = filt ;
+ PIX const *start = src + (i - filt_width) ;
+ PIX const *stop ;
+ PIX x ;
+
+ /* beginning */
+ stop = src + VL_MAX(0, i - filt_width) ;
+ x = *stop ;
+ while (start <= stop) { acc += (*g++) * x ; start++ ; }
+
+ /* middle */
+ stop = src + VL_MIN(width - 1, i + filt_width) ;
+ while (start < stop) acc += (*g++) * (*start++) ;
+
+ /* end */
+ x = *start ;
+ stop = src + (i + filt_width) ;
+ while (start <= stop) { acc += (*g++) * x ; start++ ; }
+
+ /* save */
+ *dst = acc ;
+ dst += height ;
+
+ assert (g - filt == 2 * filt_width +1) ;
+ }
+ /* next column */
+ src += width ;
+ dst -= width*height - 1 ;
+ }
+ break ;
+
+ default:
+ assert (0) ;
+ }
+}
+
+VL_EXPORT
+void
+VL_IMSMOOTH_(SFX)(PIX *dst,
+ PIX *temp,
+ PIX const *src,
+ int width, int height, double sigma)
+{
+ /* We use either a static buffer or, if the variance is very big, a
+ * new buffer. This is efficient and robust to buffer stealing. */
+ static PIX *filt = 0 ;
+ static int filt_width = -1 ;
+ static double filt_static_sigma = -1.0 ;
+ enum { filt_static_res = 1024 } ;
+ static PIX filt_static [2 * filt_static_res + 1] ;
+
+ int j ;
+
+ if (sigma < (PIX)(1e-5)) {
+ dst = memcpy(dst,src,width*height*sizeof(PIX)) ;
+ return ;
+ }
+
+ /* window width */
+ filt_width = ceil (4.0 * sigma) ;
+
+ /* setup filter only if not available from previous iteration*/
+ if (filt_static_sigma != sigma) {
+ PIX acc = 0.0 ;
+
+ if (filt_width <= filt_static_res) {
+ /* use static buffer */
+ filt = filt_static ;
+ filt_static_sigma = sigma ;
+ } else {
+ /* dynamically allocate a larger buffer */
+ filt = vl_malloc (sizeof(PIX) * (2*filt_width+1)) ;
+ }
+
+ for (j = 0 ; j < 2 * filt_width + 1 ; ++j) {
+ PIX d = (PIX)(j - filt_width) / (PIX)(sigma) ;
+ filt [j] = exp (- 0.5 * d * d) ;
+ acc += filt [j] ;
+ }
+
+ /* normalize */
+ for (j = 0 ; j < 2 * filt_width + 1 ; ++j) {
+ filt [j] /= acc ;
+ }
+ }
+
+ /* convolve */
+ VL_CONVTRANSP_(SFX) (temp, src, filt,
+ width, height, filt_width,
+ VL_CONV_CONT) ;
+ VL_CONVTRANSP_(SFX) (dst, temp, filt,
+ height, width, filt_width,
+ VL_CONV_CONT) ;
+
+ /* free buffer? */
+ if (filt_static_sigma != sigma) {
+ vl_free (filt) ;
+ }
+}
+
+/*
+ * Local Variables: *
+ * mode: C *
+ * End: *
+ */
diff --git a/vl/mathop.c b/vl/mathop.c
new file mode 100644
index 0000000..ac782e5
--- /dev/null
+++ b/vl/mathop.c
@@ -0,0 +1,30 @@
+/** @file mathop.c
+ ** @author Andrea Vedaldi
+ ** @brief Math operations - Definition
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#include "mathop.h"
+
+#define FLOAT_TYPE_FLOAT 1
+#define FLOAT_TYPE_DOUBLE 2
+
+#undef FLOAT_TYPE
+#define FLOAT_TYPE FLOAT_TYPE_DOUBLE
+
+
+#undef FLOAT_TYPE
+#define FLOAT_TYPE FLOAT_TYPE_FLOAT
+
+
+float
+vl_dist_l2_f (float * dist, int M, int NX, int NY,
+ float const* x, float const* y) ;
+
+
diff --git a/vl/mathop.h b/vl/mathop.h
new file mode 100644
index 0000000..23cf224
--- /dev/null
+++ b/vl/mathop.h
@@ -0,0 +1,450 @@
+/** @file mathop.h
+ ** @author Andrea Vedaldi
+ ** @brief Math operations
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#ifndef VL_MATHOP
+#define VL_MATHOP
+
+#include "generic.h"
+
+/** @brief Logarithm of 2 (math constant)*/
+#define VL_LOG_OF_2 0.693147180559945
+
+/** @brief Pi (math constant) */
+#define VL_PI 3.141592653589793
+
+/** @brief IEEE single precision epsilon (math constant)
+ **
+ ** <code>1.0F + VL_EPSILON_F</code> is the smallest representable
+ ** single precision number greater than @c 1.0F. Numerically,
+ ** ::VL_EPSILON_F is equal to @f$ 2^{-23} @f$.
+ **
+ **/
+#define VL_EPSILON_F 1.19209290E-07F
+
+/** @brief IEEE double precision epsilon (math constant)
+ **
+ ** <code>1.0 + VL_EPSILON_D</code> is the smallest representable
+ ** double precision number greater than @c 1.0. Numerically,
+ ** ::VL_EPSILON_D is equal to @f$ 2^{-52} @f$.
+ **/
+#define VL_EPSILON_D 2.220446049250313e-16
+
+/*
+ For the code below: An ANSI C compiler takes the two expressions,
+ LONG_VAR and CHAR_VAR, and implicitly casts them to the type of the
+ first member of the union. Refer to K&R Second Edition Page 148,
+ last paragraph.
+*/
+
+/** @brief IEEE single precision quiet NaN constant */
+static union { vl_uint32 raw ; float value ; }
+ const vl_nan_f =
+ { 0x7FC00000UL } ;
+
+/** @brief IEEE single precision infinity constant */
+static union { vl_uint32 raw ; float value ; }
+ const vl_infinity_f =
+ { 0x7F800000UL } ;
+
+/** @brief IEEE double precision quiet NaN constant */
+static union { vl_uint64 raw ; double value ; }
+ const vl_nan_d =
+#ifdef VL_COMPILER_MSC
+ { 0x7FF8000000000000ui64 } ;
+#else
+ { 0x7FF8000000000000ULL } ;
+#endif
+
+/** @brief IEEE double precision infinity constant */
+static union { vl_uint64 raw ; double value ; }
+ const vl_infinity_d =
+#ifdef VL_COMPILER_MSC
+ { 0x7FF0000000000000ui64 } ;
+#else
+ { 0x7FF0000000000000ULL } ;
+#endif
+
+/** @brief IEEE single precision NaN (not signaling) */
+#define VL_NAN_F (vl_nan_f.value)
+
+/** @brief IEEE single precision positive infinity (not signaling) */
+#define VL_INFINITY_F (vl_infinity_f.value)
+
+/** @brief IEEE double precision NaN (not signaling) */
+#define VL_NAN_D (vl_nan_d.value)
+
+/** @brief IEEE double precision positive infinity (not signaling) */
+#define VL_INFINITY_D (vl_infinity_d.value)
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast <code>mod(x, 2 * VL_PI)</code>
+ **
+ ** @param x input value.
+ **
+ ** The function is optimized for small absolute values of @a x.
+ **
+ ** The result is guaranteed not to be smaller than 0. However, due to
+ ** finite numerical precision and rounding errors, the result can be
+ ** equal to 2 * VL_PI (for instance, if @c x is a very small negative
+ ** number).
+ **
+ ** @return <code>mod(x, 2 * VL_PI)</code>
+ **/
+VL_INLINE
+float vl_mod_2pi_f (float x)
+{
+ while (x > 2 * VL_PI) x -= (float) (2 * VL_PI) ;
+ while (x < 0.0F ) x += (float) (2 * VL_PI);
+ return x ;
+}
+
+VL_INLINE
+double vl_mod_2pi_d (double x)
+{
+ while (x > 2 * VL_PI) x -= 2 * VL_PI ;
+ while (x < 0.0 ) x += 2 * VL_PI ;
+ return x ;
+}
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast <code>(int) floor(x)</code>
+ **
+ ** @param x argument.
+ ** @return @c (int) floor(x)
+ **/
+VL_INLINE
+int
+vl_floor_f (float x)
+{
+ int xi = (int) x ;
+ if (x >= 0 || (float) xi == x) return xi ;
+ else return xi - 1 ;
+}
+
+VL_INLINE
+int
+vl_floor_d (double x)
+{
+ int xi = (int) x ;
+ if (x >= 0 || (double) xi == x) return xi ;
+ else return xi - 1 ;
+}
+/* @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast @c abs(x)
+ ** @param x argument.
+ ** @return @c abs(x)
+ **/
+
+VL_INLINE
+float
+vl_abs_f (float x)
+{
+ return (x >= 0) ? x : -x ;
+}
+
+VL_INLINE
+double
+vl_abs_d (double x)
+{
+ return (x >= 0) ? x : -x ;
+}
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast @c atan2 approximation
+ ** @param y argument.
+ ** @param x argument.
+ **
+ ** The function computes a relatively rough but fast approximation of
+ ** @c atan2(y,x).
+ **
+ ** @par Algorithm
+ **
+ ** The algorithm approximates the function @f$ f(r)=atan((1-r)/(1+r))
+ ** @f$, @f$ r \in [-1,1] @f$ with a third order polynomial @f$
+ ** f(r)=c_0 + c_1 r + c_2 r^2 + c_3 r^3 @f$. To fit the polynomial
+ ** we impose the constraints
+ **
+ ** @f{eqnarray*}
+ ** f(+1) &=& c_0 + c_1 + c_2 + c_3 = atan(0) = 0,\\
+ ** f(-1) &=& c_0 - c_1 + c_2 - c_3 = atan(\infty) = \pi/2,\\
+ ** f(0) &=& c_0 = atan(1) = \pi/4.
+ ** @f}
+ **
+ ** The last degree of freedom is fixed by minimizing the @f$
+ ** l^{\infty} @f$ error, which yields
+ **
+ ** @f[
+ ** c_0=\pi/4, \quad
+ ** c_1=-0.9675, \quad
+ ** c_2=0, \quad
+ ** c_3=0.1821,
+ ** @f]
+ **
+ ** with maximum error of 0.0061 radians at 0.35 degrees.
+ **
+ ** @return Approximation of @c atan2(y,x).
+ **/
+
+VL_INLINE
+float
+vl_fast_atan2_f (float y, float x)
+{
+ float angle, r ;
+ float const c3 = 0.1821F ;
+ float const c1 = 0.9675F ;
+ float abs_y = vl_abs_f (y) + VL_EPSILON_F ;
+
+ if (x >= 0) {
+ r = (x - abs_y) / (x + abs_y) ;
+ angle = (float) (VL_PI / 4) ;
+ } else {
+ r = (x + abs_y) / (abs_y - x) ;
+ angle = (float) (3 * VL_PI / 4) ;
+ }
+ angle += (c3*r*r - c1) * r ;
+ return (y < 0) ? - angle : angle ;
+}
+
+VL_INLINE
+double
+vl_fast_atan2_d (double y, double x)
+{
+ double angle, r ;
+ double const c3 = 0.1821 ;
+ double const c1 = 0.9675 ;
+ double abs_y = vl_abs_d (y) + VL_EPSILON_D ;
+
+ if (x >= 0) {
+ r = (x - abs_y) / (x + abs_y) ;
+ angle = VL_PI / 4 ;
+ } else {
+ r = (x + abs_y) / (abs_y - x) ;
+ angle = 3 * VL_PI / 4 ;
+ }
+ angle += (c3*r*r - c1) * r ;
+ return (y < 0) ? - angle : angle ;
+}
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast @c resqrt approximation
+ ** @param x argument.
+ **
+ ** The function quickly computes an approximation of @f$ x^{-1/2}
+ ** @f$.
+ **
+ ** @par Algorithm
+ **
+ ** The goal is to compute @f$ y = x^{-1/2} @f$, which we do by
+ ** finding the solution of @f$ 0 = f(y) = y^{-2} - x@f$ by two Newton
+ ** steps. Each Newton iteration is given by
+ **
+ ** @f[
+ ** y \leftarrow
+ ** y - \frac{f(y)}{\dot f(y)} =
+ ** y + \frac{1}{2} (y-xy^3) =
+ ** \frac{y}{2} \left( 3 - xy^2 \right)
+ ** @f]
+ **
+ ** which yields a simple polynomial update rule.
+ **
+ ** The clever bit (attributed to either J. Carmack or G. Tarolli) is
+ ** the way an initial guess @f$ y \approx x^{-1/2} @f$ is chosen.
+ **
+ ** @see <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf">Inverse Sqare Root</a>.
+ **
+ ** @return Approximation of @c resqrt(x).
+ **/
+
+VL_INLINE
+float
+vl_fast_resqrt_f (float x)
+{
+ /* 32-bit version */
+ union {
+ float x ;
+ vl_int32 i ;
+ } u ;
+
+ float xhalf = (float) 0.5 * x ;
+
+ /* convert floating point value in RAW integer */
+ u.x = x ;
+
+ /* gives initial guess y0 */
+ u.i = 0x5f3759df - (u.i >> 1);
+ /*u.i = 0xdf59375f - (u.i>>1);*/
+
+ /* two Newton steps */
+ u.x = u.x * ( (float) 1.5 - xhalf*u.x*u.x) ;
+ u.x = u.x * ( (float) 1.5 - xhalf*u.x*u.x) ;
+ return u.x ;
+}
+
+VL_INLINE
+double
+vl_fast_resqrt_d (double x)
+{
+ /* 64-bit version */
+ union {
+ double x ;
+ vl_int64 i ;
+ } u ;
+
+ double xhalf = (double) 0.5 * x ;
+
+ /* convert floating point value in RAW integer */
+ u.x = x ;
+
+ /* gives initial guess y0 */
+#ifdef VL_COMPILER_MSC
+ u.i = 0x5fe6ec85e7de30dai64 - (u.i >> 1) ;
+#else
+ u.i = 0x5fe6ec85e7de30daLL - (u.i >> 1) ;
+#endif
+
+ /* two Newton steps */
+ u.x = u.x * ( (double) 1.5 - xhalf*u.x*u.x) ;
+ u.x = u.x * ( (double) 1.5 - xhalf*u.x*u.x) ;
+ return u.x ;
+}
+
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Fast @c sqrt approximation
+ ** @param x argument.
+ **
+ ** The function computes a fast approximation of @c sqrt(x).
+ **
+ ** @par Floating-point algorithm
+ **
+ ** For the floating point cases, the function uses ::vl_fast_resqrt_f
+ ** (or ::vl_fast_resqrt_d) to compute <code>x *
+ ** vl_fast_resqrt_f(x)</code>.
+ **
+ ** @par Integer algorithm
+ **
+ ** We seek for the largest integer @e y such that @f$ y^2 \leq x @f$.
+ ** Write @f$ y = w + b_k 2^k + z @f$ where the binary expansion of the
+ ** various variable is
+ **
+ ** @f[
+ ** x = \sum_{i=0}^{n-1} 2^i a_i, \qquad
+ ** w = \sum_{i=k+1}^{m-1} 2^i b_i, \qquad
+ ** z = \sum_{i=0}^{k-1} 2^i b_i.
+ ** @f]
+ **
+ ** Assume @e w known. Expanding the square and using the fact that
+ ** @f$ b_k^2=b_k @f$, we obtain the following constraint for @f$ b_k
+ ** @f$ and @e z:
+ **
+ ** @f[
+ ** x - w^2 \geq 2^k ( 2 w + 2^k ) b_k + z (z + 2wz + 2^{k+1}z b_k)
+ ** @f]
+ **
+ ** A necessary condition for @f$ b_k = 1 @f$ is that this equation is
+ ** satisfied for @f$ z = 0 @f$ (as the second term is always
+ ** non-negative). In fact, this condition is also sufficient, since
+ ** we are looking for the @e largest solution @e y.
+ **
+ ** This yields the following iterative algorithm. First, note that if
+ ** @e x is stored in @e n bits, where @e n is even, then the integer
+ ** square root does not require more than @f$ m = n / 2 @f$ bit to be
+ ** stored. Thus initially, @f$ w = 0 @f$ and @f$ k = m - 1 = n/2 - 1
+ ** @f$. Then, at each iteration the equation is tested, determining
+ ** @f$ b_{m-1}, b_{m-2}, ... @f$ in this order.
+ **
+ ** @return Approximation of @c sqrt(x).
+ **/
+
+VL_INLINE
+float
+vl_fast_sqrt_f (float x)
+{
+ return (x < 1e-8) ? 0 : x * vl_fast_resqrt_f (x) ;
+}
+
+VL_INLINE
+double
+vl_fast_sqrt_d (float x)
+{
+ return (x < 1e-8) ? 0 : x * vl_fast_resqrt_d (x) ;
+}
+
+VL_INLINE vl_uint32 vl_fast_sqrt_ui32 (vl_uint32 x) ;
+VL_INLINE vl_uint16 vl_fast_sqrt_ui16 (vl_uint16 x) ;
+VL_INLINE vl_uint8 vl_fast_sqrt_ui8 (vl_uint8 x) ;
+
+/** @} */
+
+/** @internal
+ ** @{ */
+#define VL_FAST_SQRT_UI(T, SFX) \
+ T \
+ vl_fast_sqrt_ ## SFX (T x) \
+ { \
+ T y = 0 ; /* w / 2^k */ \
+ T tmp = 0 ; \
+ int twok ; \
+ \
+ for (twok = 8 * sizeof(T) - 2 ; \
+ twok >= 0 ; twok -= 2) { \
+ y <<= 1 ; /* y = 2 * y */ \
+ tmp = (2*y + 1) << twok ; \
+ if (x >= tmp) { \
+ x -= tmp ; \
+ y += 1 ; \
+ } \
+ } \
+ return y ; \
+ }
+
+
+VL_FAST_SQRT_UI(vl_uint32, ui32)
+VL_FAST_SQRT_UI(vl_uint16, ui16)
+VL_FAST_SQRT_UI(vl_uint8, ui8 )
+/** @} */
+
+/** ------------------------------------------------------------------
+ ** @{
+ ** @brief Auto distances.
+ ** @param dist pointer to the distance matrix (out).
+ ** @param M number of dimensions
+ ** @param NX number of data points in X.
+ ** @param NY number of data points in Y.
+ ** @param X first data set.
+ ** @param Y second data set.
+ **/
+
+float
+vl_dist_l2_f (float * dist, int M, int NX, int NY,
+ float const* x, float const* y) ;
+
+float
+vl_dist_l2_d (double * dist, int M, int NX, int NY,
+ double const* x, double const* y) ;
+
+/** @} */
+
+/* VL_MATHOP */
+#endif
diff --git a/vl/mser.c b/vl/mser.c
new file mode 100644
index 0000000..4d2aea4
--- /dev/null
+++ b/vl/mser.c
@@ -0,0 +1,971 @@
+#if WIN32
+#pragma warning(disable:4018 4244)
+#endif
+
+/** @internal
+ ** @file mser.c
+ ** @brief Maximally Stable Extremal Regions (MSER) - Definition
+ ** @author Andrea Vedaldi
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+/** @file mser.h
+ **
+ ** @section mser Maximally Stable Extremal Regions Overview
+ **
+ ** Running the MSER filter usually involves the following steps:
+ **
+ ** - Initialize the MSER filter by ::vl_mser_new(). The
+ ** filter can be reused for images of the same size.
+ ** - Compute the MSERs by ::vl_mser_process().
+ ** - Optionally fit ellipsoids to the MSERs by ::vl_mser_ell_fit().
+ ** - Retrieve the results by ::vl_mser_get_regions() (and optionally ::vl_mser_get_ell()).
+ ** - Optionally retrieve filter statistics by ::vl_mser_get_stats().
+ ** - Delete the MSER filter by ::vl_mser_delete().
+ **
+ ** @subsection mser-definition MSER definition
+ **
+ ** An extremal region @f$R_l@f$ of an image is a connected component
+ ** of the level set @f$S_l = \{ x : I(x) \leq l \}@f$.
+ **
+ ** @image html mser-er.png
+ **
+ ** For each intensity @f$l@f$, one has multiple disjoint extremal
+ ** regions in the level set @f$S_l@f$. Let @f$l@f$ span a finite
+ ** number of values @f$\mathcal{L}=\{0,\dots,M-1\}@f$ (a sampling of
+ ** the image range). One obtains a family of regions @f$R_l@f$; by
+ ** connecting two regions @f$R_l@f$and @f$R_{l+1}@f$ if, and only if,
+ ** @f$R_l\subset R_{l+1}@f$, regions form a tree:
+ **
+ ** @image html mser-tree.png
+ **
+ ** The <em>maximally stable extremal regions</em> are extremal
+ ** regions which satisfy a stability criterion. Here we use a
+ ** criterion which is similar but not identical to the original
+ ** paper. This definition is somewhat simpler both to understand and
+ ** code (it also runs faster).
+ **
+ ** Let @f$B(R_l)=(R_l,R_{l+1},\dots,R_{l+\Delta})@f$ be the branch of
+ ** the tree rooted at @f$R_l@f$. We associate to the branch the
+ ** (in)stability score
+ **
+ ** @f[
+ ** v(R_l) = \frac{|R_{l+\Delta} - R_l|}{|R_l|}.
+ ** @f]
+ **
+ ** The score is low if the regions along the branch have similar area
+ ** (and thus similar shape). We aim to select maximally stable
+ ** branches; then a maximally stable region is just a representative
+ ** region selected from a maximally stable branch (for simplicity we
+ ** select @f$R_l@f$, but one could choose for example
+ ** @f$R_{l+\Delta/2}@f$).
+ **
+ ** Roughly speaking, a branch is maximally stable if it is a local
+ ** minimum of the (in)stability score. More accurately, we start by
+ ** assuming that all branches are maximally stable. Then we consider
+ ** each branch @f$B(R_{l})@f$ and its parent branch
+ ** @f$B(R_{l+1}):R_{l+1}\supset R_l@f$ (notice that, due to the
+ ** discrete nature of the calculations, they might be geometrically
+ ** identical) and we mark as unstable the less stable one, i.e.:
+ **
+ ** - if @f$v(R_l)<v(R_{l+1})@f$, mark @f$R_{l+1}@f$ as unstable;
+ ** - if @f$v(R_l)>v(R_{l+1})@f$, mark @f$R_{l}@f$ as unstable;
+ ** - otherwise, do nothing.
+ **
+ ** This criterion selects among nearby regions the ones that are more
+ ** stable. We optionally refine the selection by running (starting
+ ** from the bigger and going to the smaller regions) the following
+ ** tests:
+ **
+ ** - @f$a_- \leq |R_{l}|/|R_{\infty}| \leq a_+@f$: exclude MSERs too
+ ** small or too big (@f$|R_{\infty}|@f$ is the area of the image).
+ **
+ ** - @f$v(R_{l}) < v_+@f$: exclude MSERs too unstable.
+ **
+ ** - For any MSER @f$R_l@f$, find the parent MSER @f$R_{l'}@f$ and check
+ ** if
+ ** @f$|R_{l'} - R_l|/|R_l'| < d_+@f$: remove duplicated MSERs.
+ **
+ ** <table>
+ ** <tr>
+ ** <td>parameter</td>
+ ** <td>alt. name</td>
+ ** <td>standard value</td>
+ ** <td>set by</td>
+ ** </tr>
+ ** <tr>
+ ** <td>@f$\Delta@f$</td>
+ ** <td>@c delta</td>
+ ** <td>5</td>
+ ** <td>::vl_mser_set_delta()</td>
+ ** </tr>
+ ** <tr>
+ ** <td>@f$a_+@f$</td>
+ ** <td>@c max_area</td>
+ ** <td>0.75</td>
+ ** <td>::vl_mser_set_max_area()</td>
+ ** </tr>
+ ** <tr>
+ ** <td>@f$a_-@f$</td>
+ ** <td>@c min_area</td>
+ ** <td>3.0/@f$|R_\infty|@f$</td>
+ ** <td>::vl_mser_set_min_area()</td>
+ ** </tr>
+ ** <tr>
+ ** <td>@f$v_+@f$</td>
+ ** <td>@c max_var</td>
+ ** <td>0.25</td>
+ ** <td>::vl_mser_set_max_variation()</td>
+ ** </tr>
+ ** <tr>
+ ** <td>@f$d_+@f$</td>
+ ** <td>@c min_diversity</td>
+ ** <td>0.2</td>
+ ** <td>::vl_mser_set_min_diversity()</td>
+ ** </tr>
+ ** </table>
+ **
+ ** @subsection mser-vol Volumetric images
+ **
+ ** The code supports images of arbitrary dimension. For instance, it
+ ** is possible to find the MSER regions of volumetric images or time
+ ** sequences. See ::vl_mser_new() for further details.s
+ **
+ ** @subsection mser-ell Ellipsoids
+ **
+ ** Usually extremal regions are returned as a set of ellipsoids
+ ** fitted to the actual regions (which have arbitrary shape). The fit
+ ** is done by calculating the mean and variance of the pixels
+ ** composing the region:
+ ** @f[
+ ** \mu_l = \frac{1}{|R_l|}\sum_{x\in R_l}x,
+ ** \qquad
+ ** \Sigma_l = \frac{1}{|R_l|}\sum_{x\in R_l} (x-\mu_l)^\top(x-\mu_l)
+ ** @f]
+ ** Ellipsoids are fitted by ::vl_mser_ell_fit(). Notice that for a
+ ** <em>n</em> dimensional image, the mean has <em>n</em> components
+ ** and the variance has <em>n(n+1)/2</em> independent components. The
+ ** total number of components is obtained by ::vl_mser_get_ell_dof()
+ ** and the total number of fitted ellipsoids by
+ ** ::vl_mser_get_ell_num(). A matrix with an ellipsoid per column is
+ ** returned by ::vl_mser_get_ell(). The column is the stacking of the
+ ** mean and of the independent components of the variance, in the
+ ** order <em>(1,1),(1,2),..,(1,n), (2,2),(2,3)...</em>. In the
+ ** calculations, the pixel coordinate @f$x=(x_1,...,x_n)@f$ use the
+ ** standard index order and ranges.
+ **
+ ** @subsection mser-algo Algorithm
+ **
+ ** The algorithm is quite efficient. While some details may be
+ ** tricky, the overall idea is easy to grasp.
+ **
+ ** - Pixels are sorted by increasing intensity.
+ ** - Pixels are added to a forest by increasing intensity. The forest has the
+ ** following properties:
+ ** - All the descendent of a certain pixels are subset of an extremal region.
+ ** - All the extremal regions are the descendants of some pixels.
+ ** - Extremal regions are extracted from the region tree and the extremal regions tree is
+ ** calculated.
+ ** - Stable regions are marked.
+ ** - Duplicates and other bad regions are removed.
+ **
+ ** @remark The extremal region tree which is calculated is a subset
+ ** of the actual extremal region tree. In particular, it does not
+ ** contain redundant entries extremal regions that coincide as
+ ** sets. So, for example, in the calculated extremal region tree, the
+ ** parent @f$R_q@f$ of an extremal region @f$R_{l}@f$ may or may
+ ** <em>not</em> correspond to @f$R_{l+1}@f$, depending whether
+ ** @f$q\leq l+1@f$ or not. These subtleties are important when
+ ** calculating the stability tests.
+ **/
+
+#include "mser.h"
+
+#include<stdlib.h>
+#include<string.h>
+#include<assert.h>
+
+/** -------------------------------------------------------------------
+ ** @brief Advance N-dimensional subscript
+ **
+ ** The function increments by one the subscript @a subs indexing an
+ ** array the @a ndims dimensions @a dims.
+ **
+ ** @param ndims number of dimensions.
+ ** @param dims dimensions.
+ ** @param subs subscript to advance.
+ **/
+
+VL_INLINE void
+adv(int ndims, int const *dims, int *subs)
+{
+ int d = 0 ;
+ while(d < ndims) {
+ if( ++subs[d] < dims[d] ) return ;
+ subs[d++] = 0 ;
+ }
+}
+
+/** -------------------------------------------------------------------
+ ** @brief Climb the region forest to reach aa root
+ **
+ ** The function climbs the regions forest @a r starting from the node
+ ** @a idx to the corresponding root.
+ **
+ ** To speed-up the operation, the function uses the
+ ** VlMserReg::shortcut field to quickly jump to the root. After the
+ ** root is reached, all the used shortcut are updated.
+ **
+ ** @param r regions' forest.
+ ** @param idx stating node.
+ ** @return index of the reached root.
+ **/
+
+VL_INLINE vl_uint
+climb (VlMserReg* r, vl_uint idx)
+{
+
+ vl_uint prev_idx = idx ;
+ vl_uint next_idx ;
+ vl_uint root_idx ;
+
+ /* move towards root to find it */
+ while (1) {
+
+ /* next jump to the root */
+ next_idx = r [idx] .shortcut ;
+
+ /* recycle shortcut to remember how we came here */
+ r [idx] .shortcut = prev_idx ;
+
+ /* stop if the root is found */
+ if( next_idx == idx ) break ;
+
+ /* next guy */
+ prev_idx = idx ;
+ idx = next_idx ;
+ }
+
+ root_idx = idx ;
+
+ /* move backward to update shortcuts */
+ while (1) {
+
+ /* get previously visited one */
+ prev_idx = r [idx] .shortcut ;
+
+ /* update shortcut to point to the new root */
+ r [idx] .shortcut = root_idx ;
+
+ /* stop if the first visited node is reached */
+ if( prev_idx == idx ) break ;
+
+ /* next guy */
+ idx = prev_idx ;
+ }
+
+ return root_idx ;
+}
+
+/** -------------------------------------------------------------------
+ ** @brief Create a new MSER filter
+ **
+ ** Initializes a new MSER filter for images of the specified
+ ** dimensions. Images are @a ndims -dimensional arrays of dimensions
+ ** @a dims.
+ **
+ ** @param ndims number of dimensions.
+ ** @param dims dimensions.
+ **/
+VL_EXPORT
+VlMserFilt*
+vl_mser_new (int ndims, int const* dims)
+{
+ VlMserFilt* f ;
+ int *strides, k ;
+
+ f = vl_calloc (sizeof(VlMserFilt), 1) ;
+
+ f-> ndims = ndims ;
+ f-> dims = vl_malloc (sizeof(int) * ndims) ;
+ f-> subs = vl_malloc (sizeof(int) * ndims) ;
+ f-> dsubs = vl_malloc (sizeof(int) * ndims) ;
+ f-> strides = vl_malloc (sizeof(int) * ndims) ;
+
+ /* shortcuts */
+ strides = f-> strides ;
+
+ /* copy dims to f->dims */
+ for(k = 0 ; k < ndims ; ++k) {
+ f-> dims [k] = dims [k] ;
+ }
+
+ /* compute strides to move into the N-dimensional image array */
+ strides [0] = 1 ;
+ for(k = 1 ; k < ndims ; ++k) {
+ strides [k] = strides [k-1] * dims [k-1] ;
+ }
+
+ /* total number of pixels */
+ f-> nel = strides [ndims-1] * dims [ndims-1] ;
+
+ /* dof of ellipsoids */
+ f-> dof = ndims * (ndims + 1) / 2 + ndims ;
+
+ /* more buffers */
+ f-> perm = vl_malloc (sizeof(vl_uint) * f-> nel) ;
+ f-> joins = vl_malloc (sizeof(vl_uint) * f-> nel) ;
+ f-> r = vl_malloc (sizeof(VlMserReg) * f-> nel) ;
+
+ f-> er = 0 ;
+ f-> rer = 0 ;
+ f-> mer = 0 ;
+ f-> rmer = 0 ;
+ f-> ell = 0 ;
+ f-> rell = 0 ;
+
+ /* other parameters */
+ f-> delta = 5 ;
+ f-> max_area = 0.75 ;
+ f-> min_area = 3.0 / f-> nel ;
+ f-> max_variation = 0.25 ;
+ f-> min_diversity = 0.2 ;
+
+ return f ;
+}
+
+/** -------------------------------------------------------------------
+ ** @brief Delete MSER filter
+ **
+ ** The function releases the MSER filter @a f and all its resources.
+ **
+ ** @param f MSER filter to be deleted.
+ **/
+VL_EXPORT
+void
+vl_mser_delete (VlMserFilt* f)
+{
+ if(f) {
+ if(f-> acc ) vl_free( f-> acc ) ;
+ if(f-> ell ) vl_free( f-> ell ) ;
+
+ if(f-> er ) vl_free( f-> er ) ;
+ if(f-> r ) vl_free( f-> r ) ;
+ if(f-> joins ) vl_free( f-> joins ) ;
+ if(f-> perm ) vl_free( f-> perm ) ;
+
+ if(f-> strides) vl_free( f-> strides) ;
+ if(f-> dsubs ) vl_free( f-> dsubs ) ;
+ if(f-> subs ) vl_free( f-> subs ) ;
+ if(f-> dims ) vl_free( f-> dims ) ;
+
+ if(f-> mer ) vl_free( f-> mer ) ;
+ vl_free (f) ;
+ }
+}
+
+
+/** -------------------------------------------------------------------
+ ** @brief Process image
+ **
+ ** The functions calculates the Maximally Stable Extremal Regions
+ ** (MSERs) of image @a im using the MSER filter @a f.
+ **
+ ** The filter @a f must have been initialized to be compatible with
+ ** the dimensions of @a im.
+ **
+ ** @param f MSER filter.
+ ** @param im image data.
+ **/
+VL_EXPORT
+void
+vl_mser_process (VlMserFilt* f, vl_mser_pix const* im)
+{
+ /* shortcuts */
+ vl_uint nel = f-> nel ;
+ vl_uint *perm = f-> perm ;
+ vl_uint *joins = f-> joins ;
+ int ndims = f-> ndims ;
+ int *dims = f-> dims ;
+ int *subs = f-> subs ;
+ int *dsubs = f-> dsubs ;
+ int *strides = f-> strides ;
+ VlMserReg *r = f-> r ;
+ VlMserExtrReg *er = f-> er ;
+ vl_uint *mer = f-> mer ;
+ int delta = f-> delta ;
+
+ int njoins = 0 ;
+ int ner = 0 ;
+ int nmer = 0 ;
+ int nbig = 0 ;
+ int nsmall = 0 ;
+ int nbad = 0 ;
+ int ndup = 0 ;
+
+ int i, j, k ;
+
+ /* delete any previosuly computed ellipsoid */
+ f-> nell = 0 ;
+
+ /* -----------------------------------------------------------------
+ * Sort pixels by intensity
+ * -------------------------------------------------------------- */
+
+ {
+ vl_uint buckets [ VL_MSER_PIX_MAXVAL ] ;
+
+ /* clear buckets */
+ memset (buckets, 0, sizeof(vl_uint) * VL_MSER_PIX_MAXVAL ) ;
+
+ /* compute bucket size (how many pixels for each intensity
+ value) */
+ for(i = 0 ; i < nel ; ++i) {
+ vl_mser_pix v = im [i] ;
+ ++ buckets [v] ;
+ }
+
+ /* cumulatively add bucket sizes */
+ for(i = 1 ; i < VL_MSER_PIX_MAXVAL ; ++i) {
+ buckets [i] += buckets [i-1] ;
+ }
+
+ /* empty buckets computing pixel ordering */
+ for(i = nel ; i >= 1 ; ) {
+ vl_mser_pix v = im [ --i ] ;
+ vl_uint j = -- buckets [v] ;
+ perm [j] = i ;
+ }
+ }
+
+ /* initialize the forest with all void nodes */
+ for(i = 0 ; i < nel ; ++i) {
+ r [i] .parent = VL_MSER_VOID_NODE ;
+ }
+
+ /* -----------------------------------------------------------------
+ * Compute regions and count extremal regions
+ * -------------------------------------------------------------- */
+ /*
+ In the following:
+
+ idx : index of the current pixel
+ val : intensity of the current pixel
+ r_idx : index of the root of the current pixel
+ n_idx : index of the neighbors of the current pixel
+ nr_idx : index of the root of the neighbor of the current pixel
+
+ */
+
+ /* process each pixel by increasing intensity */
+ for(i = 0 ; i < nel ; ++i) {
+
+ /* pop next node xi */
+ vl_uint idx = perm [i] ;
+ vl_mser_pix val = im [idx] ;
+ vl_uint r_idx ;
+
+ /* add the pixel to the forest as a root for now */
+ r [idx] .parent = idx ;
+ r [idx] .shortcut = idx ;
+ r [idx] .area = 1 ;
+ r [idx] .height = 1 ;
+
+ r_idx = idx ;
+
+ /* convert the index IDX into the subscript SUBS; also initialize
+ DSUBS to (-1,-1,...,-1) */
+ {
+ vl_uint temp = idx ;
+ for(k = ndims - 1 ; k >= 0 ; --k) {
+ dsubs [k] = -1 ;
+ subs [k] = temp / strides [k] ;
+ temp = temp % strides [k] ;
+ }
+ }
+
+ /* examine the neighbors of the current pixel */
+ while (1) {
+ vl_uint n_idx = 0 ;
+ vl_bool good = 1 ;
+
+ /*
+ Compute the neighbor subscript as NSUBS+SUB, the
+ corresponding neighbor index NINDEX and check that the
+ neighbor is within the image domain.
+ */
+ for(k = 0 ; k < ndims && good ; ++k) {
+ int temp = dsubs [k] + subs [k] ;
+ good &= (0 <= temp) && (temp < dims [k]) ;
+ n_idx += temp * strides [k] ;
+ }
+
+ /*
+ The neighbor should be processed if the following conditions
+ are met:
+
+ 1. The neighbor is within image boundaries.
+
+ 2. The neighbor is indeed different from the current node
+ (the opposite happens when DSUB=(0,0,...,0)).
+
+ 3. The neighbor is already in the forest, meaning that it has
+ already been processed.
+ */
+ if (good &&
+ n_idx != idx &&
+ r [n_idx] .parent != VL_MSER_VOID_NODE ) {
+
+ vl_mser_pix nr_val = 0 ;
+ vl_uint nr_idx = 0 ;
+ int hgt = r [ r_idx] .height ;
+ int n_hgt = r [nr_idx] .height ;
+
+ /*
+ Now we join the two subtrees rooted at
+
+ R_IDX = ROOT( IDX)
+ NR_IDX = ROOT(N_IDX).
+
+ Note that R_IDX = ROOT(IDX) might change as we process more
+ neighbors, so we need keep updating it.
+ */
+
+ r_idx = climb(r, idx) ;
+ nr_idx = climb(r, n_idx) ;
+
+ /*
+ At this point we have three possibilities:
+
+ (A) ROOT(IDX) == ROOT(NR_IDX). In this case the two trees
+ have already been joined and we do not do anything.
+
+ (B) I(ROOT(IDX)) == I(ROOT(NR_IDX)). In this case the pixel
+ IDX is extending an extremal region with the same
+ intensity value. Since ROOT(NR_IDX) will NOT be an
+ extremal region of the full image, ROOT(IDX) can be
+ safely added as children of ROOT(NR_IDX) if this
+ reduces the height according to the union rank
+ heuristic.
+
+ (C) I(ROOT(IDX)) > I(ROOT(NR_IDX)). In this case the pixel
+ IDX is starting a new extremal region. Thus ROOT(NR_IDX)
+ WILL be an extremal region of the final image and the
+ only possibility is to add ROOT(NR_IDX) as children of
+ ROOT(IDX), which becomes parent.
+ */
+
+ if( r_idx != nr_idx ) { /* skip if (A) */
+
+ nr_val = im [nr_idx] ;
+
+ if( nr_val == val && hgt < n_hgt ) {
+
+ /* ROOT(IDX) becomes the child */
+ r [r_idx] .parent = nr_idx ;
+ r [r_idx] .shortcut = nr_idx ;
+ r [nr_idx] .area += r [r_idx] .area ;
+ r [nr_idx] .height = VL_MAX(n_hgt, hgt+1) ;
+
+ joins [njoins++] = r_idx ;
+
+ } else {
+
+ /* cases ROOT(IDX) becomes the parent */
+ r [nr_idx] .parent = r_idx ;
+ r [nr_idx] .shortcut = r_idx ;
+ r [r_idx] .area += r [nr_idx] .area ;
+ r [r_idx] .height = VL_MAX(hgt, n_hgt + 1) ;
+
+ joins [njoins++] = nr_idx ;
+
+ /* count if extremal */
+ if (nr_val != val) ++ ner ;
+
+ } /* check b vs c */
+ } /* check a vs b or c */
+ } /* neighbor done */
+
+ /* move to next neighbor */
+ k = 0 ;
+ while(++ dsubs [k] > 1) {
+ dsubs [k++] = -1 ;
+ if(k == ndims) goto done_all_neighbors ;
+ }
+ } /* next neighbor */
+ done_all_neighbors : ;
+ } /* next pixel */
+
+ /* the last root is extremal too */
+ ++ ner ;
+
+ /* save back */
+ f-> njoins = njoins ;
+
+ f-> stats. num_extremal = ner ;
+
+ /* -----------------------------------------------------------------
+ * Extract extremal regions
+ * -------------------------------------------------------------- */
+
+ /*
+ Extremal regions are extracted and stored into the array ER. The
+ structure R is also updated so that .SHORTCUT indexes the
+ corresponding extremal region if any (otherwise it is set to
+ VOID).
+ */
+
+ /* make room */
+ if (f-> rer < ner) {
+ if (er) vl_free (er) ;
+ f->er = er = vl_malloc (sizeof(VlMserExtrReg) * ner) ;
+ f->rer = ner ;
+ } ;
+
+ /* save back */
+ f-> nmer = ner ;
+
+ /* count again */
+ ner = 0 ;
+
+ /* scan all regions Xi */
+ for(i = 0 ; i < nel ; ++i) {
+
+ /* pop next node xi */
+ vl_uint idx = perm [i] ;
+
+ vl_mser_pix val = im [idx] ;
+ vl_uint p_idx = r [idx] .parent ;
+ vl_mser_pix p_val = im [p_idx] ;
+
+ /* is extremal ? */
+ vl_bool is_extr = (p_val > val) || idx == p_idx ;
+
+ if( is_extr ) {
+
+ /* if so, add it */
+ er [ner] .index = idx ;
+ er [ner] .parent = ner ;
+ er [ner] .value = im [idx] ;
+ er [ner] .area = r [idx] .area ;
+
+ /* link this region to this extremal region */
+ r [idx] .shortcut = ner ;
+
+ /* increase count */
+ ++ ner ;
+ } else {
+ /* link this region to void */
+ r [idx] .shortcut = VL_MSER_VOID_NODE ;
+ }
+ }
+
+ /* -----------------------------------------------------------------
+ * Link extremal regions in a tree
+ * -------------------------------------------------------------- */
+
+ for(i = 0 ; i < ner ; ++i) {
+
+ vl_uint idx = er [i] .index ;
+
+ do {
+ idx = r[idx] .parent ;
+ } while (r[idx] .shortcut == VL_MSER_VOID_NODE) ;
+
+ er[i] .parent = r[idx] .shortcut ;
+ er[i] .shortcut = i ;
+ }
+
+ /* -----------------------------------------------------------------
+ * Compute variability of +DELTA branches
+ * -------------------------------------------------------------- */
+ /* For each extremal region Xi of value VAL we look for the biggest
+ * parent that has value not greater than VAL+DELTA. This is dubbed
+ * `top parent'. */
+
+ for(i = 0 ; i < ner ; ++i) {
+
+ /* Xj is the current region the region and Xj are the parents */
+ int top_val = er [i] .value + delta ;
+ int top = er [i] .shortcut ;
+
+ /* examine all parents */
+ while (1) {
+ int next = er [top] .parent ;
+ int next_val = er [next] .value ;
+
+ /* Break if:
+ * - there is no node above the top or
+ * - the next node is above the top value.
+ */
+ if (next == top || next_val > top_val) break ;
+
+ /* so next could be the top */
+ top = next ;
+ }
+
+ /* calculate branch variation */
+ {
+ int area = er [i ] .area ;
+ int area_top = er [top] .area ;
+ er [i] .variation = (float) (area_top - area) / area ;
+ er [i] .max_stable = 1 ;
+ }
+
+ /* Optimization: since extremal regions are processed by
+ * increasing intensity, all next extremal regions being processed
+ * have value at least equal to the one of Xi. If any of them has
+ * parent the parent of Xi (this comprises the parent itself), we
+ * can safely skip most intermediate node along the branch and
+ * skip directly to the top to start our search. */
+ {
+ int parent = er [i] .parent ;
+ int curr = er [parent] .shortcut ;
+ er [parent] .shortcut = VL_MAX (top, curr) ;
+ }
+ }
+
+ /* -----------------------------------------------------------------
+ * Select maximally stable branches
+ * -------------------------------------------------------------- */
+
+ nmer = ner ;
+ for(i = 0 ; i < ner ; ++i) {
+ vl_uint parent = er [i ] .parent ;
+ vl_mser_pix val = er [i ] .value ;
+ float var = er [i ] .variation ;
+ vl_mser_pix p_val = er [parent] .value ;
+ float p_var = er [parent] .variation ;
+ vl_uint loser ;
+
+ /*
+ Notice that R_parent = R_{l+1} only if p_val = val + 1. If not,
+ this and the parent region coincide and there is nothing to do.
+ */
+ if(p_val > val + 1) continue ;
+
+ /* decide which one to keep and put that in loser */
+ if(var < p_var) loser = parent ; else loser = i ;
+
+ /* make loser NON maximally stable */
+ if(er [loser] .max_stable) {
+ -- nmer ;
+ er [loser] .max_stable = 0 ;
+ }
+ }
+
+ f-> stats. num_unstable = ner - nmer ;
+
+ /* -----------------------------------------------------------------
+ * Further filtering
+ * -------------------------------------------------------------- */
+ /* It is critical for correct duplicate detection to remove regions
+ * from the bottom (smallest one first). */
+ {
+ float max_area = (float) f-> max_area * nel ;
+ float min_area = (float) f-> min_area * nel ;
+ float max_var = (float) f-> max_variation ;
+ float min_div = (float) f-> min_diversity ;
+
+ /* scan all extremal regions (intensity value order) */
+ for(i = ner-1 ; i >= 0L ; --i) {
+
+ /* process only maximally stable extremal regions */
+ if (! er [i] .max_stable) continue ;
+
+ if (er [i] .variation >= max_var ) { ++ nbad ; goto remove ; }
+ if (er [i] .area > max_area) { ++ nbig ; goto remove ; }
+ if (er [i] .area < min_area) { ++ nsmall ; goto remove ; }
+
+ /*
+ * Remove duplicates
+ */
+ if (min_div < 1.0) {
+ vl_uint parent = er [i] .parent ;
+ int area, p_area ;
+ float div ;
+
+ /* check all but the root mser */
+ if(parent != i) {
+
+ /* search for the maximally stable parent region */
+ while(! er [parent] .max_stable) {
+ vl_uint next = er [parent] .parent ;
+ if(next == parent) break ;
+ parent = next ;
+ }
+
+ /* Compare with the parent region; if the current and parent
+ * regions are too similar, keep only the parent. */
+ area = er [i] .area ;
+ p_area = er [parent] .area ;
+ div = (float) (p_area - area) / (float) p_area ;
+
+ if (div < min_div) { ++ ndup ; goto remove ; }
+ } /* remove dups end */
+
+ }
+ continue ;
+ remove :
+ er [i] .max_stable = 0 ;
+ -- nmer ;
+ } /* check next region */
+
+ f-> stats .num_abs_unstable = nbad ;
+ f-> stats .num_too_big = nbig ;
+ f-> stats .num_too_small = nsmall ;
+ f-> stats .num_duplicates = ndup ;
+ }
+ /* -----------------------------------------------------------------
+ * Save the result
+ * -------------------------------------------------------------- */
+
+ /* make room */
+ if (f-> rmer < nmer) {
+ if (mer) vl_free (mer) ;
+ f->mer = mer = vl_malloc( sizeof(vl_uint) * nmer) ;
+ f->rmer = nmer ;
+ }
+
+ /* save back */
+ f-> nmer = nmer ;
+
+ j = 0 ;
+ for (i = 0 ; i < ner ; ++i) {
+ if (er [i] .max_stable) mer [j++] = er [i] .index ;
+ }
+}
+
+/** -------------------------------------------------------------------
+ ** @brief Fit ellipsoids
+ **
+ ** @param f MSER filter.
+ **
+ ** @sa @ref mser-ell
+ **/
+
+VL_EXPORT
+void
+vl_mser_ell_fit (VlMserFilt* f)
+{
+ /* shortcuts */
+ int nel = f-> nel ;
+ int dof = f-> dof ;
+ int *dims = f-> dims ;
+ int ndims = f-> ndims ;
+ int *subs = f-> subs ;
+ int njoins = f-> njoins ;
+ vl_uint *joins = f-> joins ;
+ VlMserReg *r = f-> r ;
+ vl_uint *mer = f-> mer ;
+ int nmer = f-> nmer ;
+ vl_mser_acc *acc = f-> acc ;
+ vl_mser_acc *ell = f-> ell ;
+
+ int d, index, i, j ;
+
+ /* already fit ? */
+ if (f->nell == f->nmer) return ;
+
+ /* make room */
+ if (f->rell < f->nmer) {
+ if (f->ell) vl_free (f->ell) ;
+ f->ell = vl_malloc (sizeof(float) * f->nmer * f->dof) ;
+ f->rell = f-> nmer ;
+ }
+
+ if (f->acc == 0) {
+ f->acc = vl_malloc (sizeof(float) * f->nel) ;
+ }
+
+ acc = f-> acc ;
+ ell = f-> ell ;
+
+ /* -----------------------------------------------------------------
+ * Integrate moments
+ * -------------------------------------------------------------- */
+
+ /* for each dof */
+ for(d = 0 ; d < f->dof ; ++d) {
+
+ /* start from the upper-left pixel (0,0,...,0) */
+ memset (subs, 0, sizeof(int) * ndims) ;
+
+ /* step 1: fill acc pretending that each region has only one pixel */
+ if(d < ndims) {
+ /* 1-order ................................................... */
+
+ for(index = 0 ; index < nel ; ++ index) {
+ acc [index] = subs [d] ;
+ adv(ndims, dims, subs) ;
+ }
+ }
+ else {
+ /* 2-order ................................................... */
+
+ /* map the dof d to a second order moment E[x_i x_j] */
+ i = d - ndims ;
+ j = 0 ;
+ while(i > j) {
+ i -= j + 1 ;
+ j ++ ;
+ }
+ /* initialize acc with x_i * x_j */
+ for(index = 0 ; index < nel ; ++ index){
+ acc [index] = subs [i] * subs [j] ;
+ adv(ndims, dims, subs) ;
+ }
+ }
+
+ /* step 2: integrate */
+ for(i = 0 ; i < njoins ; ++i) {
+ vl_uint index = joins [i] ;
+ vl_uint parent = r [index] .parent ;
+ acc [parent] += acc [index] ;
+ }
+
+ /* step 3: save back to ellpises */
+ for(i = 0 ; i < nmer ; ++i) {
+ vl_uint idx = mer [i] ;
+ ell [d + dof*i] = acc [idx] ;
+ }
+
+ } /* next dof */
+
+ /* -----------------------------------------------------------------
+ * Compute central moments
+ * -------------------------------------------------------------- */
+
+ for(index = 0 ; index < nmer ; ++index) {
+ float *pt = ell + index * dof ;
+ vl_uint idx = mer [index] ;
+ float area = r [idx] .area ;
+
+ for(d = 0 ; d < dof ; ++d) {
+
+ pt [d] /= area ;
+
+ if(d >= ndims) {
+ /* remove squared mean from moment to get variance */
+ i = d - ndims ;
+ j = 0 ;
+ while(i > j) {
+ i -= j + 1 ;
+ j ++ ;
+ }
+ pt [d] -= pt [i] * pt [j] ;
+ }
+
+ }
+ }
+
+ /* save back */
+ f-> nell = nmer ;
+}
+
diff --git a/vl/mser.h b/vl/mser.h
new file mode 100644
index 0000000..80327e1
--- /dev/null
+++ b/vl/mser.h
@@ -0,0 +1,422 @@
+/** @file mser.h
+ ** @brief Maximally Stable Extremal Regions (MSER)
+ ** @author Andrea Vedaldi
+ **
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#ifndef VL_MSER
+#define VL_MSER
+
+#include "generic.h"
+
+/** @brief MSER image data type
+ **
+ ** This is the data type of the image pixels. It has to be an
+ ** integer.
+ **/
+typedef vl_uint8 vl_mser_pix ;
+
+/** @brief Maximum value
+ **
+ ** Maximum value of the integer type ::vl_mser_pix.
+ **/
+#define VL_MSER_PIX_MAXVAL 256
+
+
+/** @brief MSER Filter
+ **
+ ** The MSER filter computes the Maximally Stable Extremal Regions of
+ ** an image.
+ **
+ ** This structure is @ref main-glossary "opaque".
+ **
+ ** @sa @ref mser
+ **/
+typedef struct _VlMserFilt VlMserFilt ;
+
+/** @brief MSER filter statistics */
+typedef struct _VlMserStats VlMserStats ;
+
+/** @brief MSER filter statistics definition */
+struct _VlMserStats
+{
+ int num_extremal ; /**< number of extremal regions */
+ int num_unstable ; /**< number of unstable extremal regions */
+ int num_abs_unstable ; /**< number of regions that failed the absolute stability test */
+ int num_too_big ; /**< number of regions that failed the maximum size test */
+ int num_too_small ; /**< number of regions that failed the minimum size test */
+ int num_duplicates ; /**< number of regions that failed the duplicate test */
+} ;
+
+/** @name Construction and Destruction
+ ** @{
+ **/
+VL_EXPORT VlMserFilt* vl_mser_new (int ndims, int const* dims) ;
+VL_EXPORT void vl_mser_delete (VlMserFilt *f) ;
+/** @} */
+
+/** @name Processing
+ ** @{
+ **/
+VL_EXPORT void vl_mser_process (VlMserFilt *f,
+ vl_mser_pix const *im) ;
+VL_EXPORT void vl_mser_ell_fit (VlMserFilt *f) ;
+/** @} */
+
+/** @name Retrieving data
+ ** @{
+ **/
+VL_INLINE vl_uint vl_mser_get_regions_num (VlMserFilt const *f) ;
+VL_INLINE vl_uint const* vl_mser_get_regions (VlMserFilt const *f) ;
+VL_INLINE float const* vl_mser_get_ell (VlMserFilt const *f) ;
+VL_INLINE vl_uint vl_mser_get_ell_num (VlMserFilt const *f) ;
+VL_INLINE vl_uint vl_mser_get_ell_dof (VlMserFilt const *f) ;
+VL_INLINE VlMserStats const* vl_mser_get_stats (VlMserFilt const *f) ;
+/** @} */
+
+/** @name Retrieving parameters
+ ** @{
+ **/
+VL_INLINE vl_mser_pix vl_mser_get_delta (VlMserFilt const *f) ;
+VL_INLINE double vl_mser_get_min_area (VlMserFilt const *f) ;
+VL_INLINE double vl_mser_get_max_area (VlMserFilt const *f) ;
+VL_INLINE double vl_mser_get_max_variation (VlMserFilt const *f) ;
+VL_INLINE double vl_mser_get_min_diversity (VlMserFilt const *f) ;
+/** @} */
+
+/** @name Setting parameters
+ ** @{
+ **/
+VL_INLINE void vl_mser_set_delta (VlMserFilt *f, vl_mser_pix x) ;
+VL_INLINE void vl_mser_set_min_area (VlMserFilt *f, double x) ;
+VL_INLINE void vl_mser_set_max_area (VlMserFilt *f, double x) ;
+VL_INLINE void vl_mser_set_max_variation (VlMserFilt *f, double x) ;
+VL_INLINE void vl_mser_set_min_diversity (VlMserFilt *f, double x) ;
+/** @} */
+
+/* ====================================================================
+ * INLINE DEFINITIONS
+ * ================================================================== */
+
+/** @internal @brief MSER accumulator data type
+ **
+ ** This is a large integer type. It should be large enough to contain
+ ** a number equal to the area (volume) of the image by the image
+ ** width by the image height (for instance, if the image is a square
+ ** of side 256, the maximum value is 256 x 256 x 256).
+ **/
+typedef float vl_mser_acc ;
+
+/** @internal @brief Basic region flag: null region */
+#ifdef VL_COMPILER_MSC
+#define VL_MSER_VOID_NODE ((1ui64<<32) - 1)
+#else
+#define VL_MSER_VOID_NODE ((1ULL<<32) - 1)
+#endif
+
+/* ----------------------------------------------------------------- */
+/** @internal @brief MSER: basic region (declaration)
+ **
+ ** Extremal regions and maximally stable extremal regions are
+ ** instances of image regions.
+ **
+ ** There is an image region for each pixel of the image. Each region
+ ** is represented by an instance of this structure. Regions are
+ ** stored into an array in pixel order.
+ **
+ ** Regions are arranged into a forest. VlMserReg::parent points to
+ ** the parent node, or to the node itself if the node is a root.
+ ** VlMserReg::parent is the index of the node in the node array
+ ** (which therefore is also the index of the corresponding
+ ** pixel). VlMserReg::height is the distance of the fartest leaf. If
+ ** the node itself is a leaf, then VlMserReg::height is zero.
+ **
+ ** VlMserReg::area is the area of the image region corresponding to
+ ** this node.
+ **
+ ** VlMserReg::region is the extremal region identifier. Not all
+ ** regions are extremal regions however; if the region is NOT
+ ** extremal, this field is set to ....
+ **/
+struct _VlMserReg
+{
+ vl_uint parent ; /**< points to the parent region. */
+ vl_uint shortcut ; /**< points to a region closer to a root. */
+ vl_uint height ; /**< region height in the forest. */
+ vl_uint area ; /**< area of the region. */
+} ;
+
+/** @internal @brief MSER: basic region */
+typedef struct _VlMserReg VlMserReg ;
+
+/* ----------------------------------------------------------------- */
+/** @internal @brief MSER: extremal region (declaration)
+ **
+ ** Extremal regions (ER) are extracted from the region forest. Each
+ ** region is represented by an instance of this structure. The
+ ** structures are stored into an array, in arbitrary order.
+ **
+ ** ER are arranged into a tree. @a parent points to the parent ER, or
+ ** to itself if the ER is the root.
+ **
+ ** An instance of the structure represents the extremal region of the
+ ** level set of intensity VlMserExtrReg::value and containing the
+ ** pixel VlMserExtReg::index.
+ **
+ ** VlMserExtrReg::area is the area of the extremal region and
+ ** VlMserExtrReg::area_top is the area of the extremal region
+ ** containing this region in the level set of intensity
+ ** VlMserExtrReg::area + @c delta.
+ **
+ ** VlMserExtrReg::variation is the relative area variation @c
+ ** (area_top-area)/area.
+ **
+ ** VlMserExtrReg::max_stable is a flag signaling whether this extremal
+ ** region is also maximally stable.
+ **/
+struct _VlMserExtrReg
+{
+ int parent ; /**< index of the parent region */
+ int index ; /**< index of pivot pixel */
+ vl_mser_pix value ; /**< value of pivot pixel */
+ vl_uint shortcut ; /**< shortcut used when building a tree */
+ vl_uint area ; /**< area of the region */
+ float variation ; /**< rel. area variation */
+ vl_uint max_stable ; /**< max stable number (=0 if not maxstable) */
+} ;
+
+/** @internal @brief MSER: extremal region */
+typedef struct _VlMserExtrReg VlMserExtrReg ;
+
+/* ----------------------------------------------------------------- */
+/** @internal @brief MSER filter
+ ** @see @ref mser
+ **/
+struct _VlMserFilt
+{
+
+ /** @name Image data and meta data @internal */
+ /*@{*/
+ int ndims ; /**< number of dimensions */
+ int *dims ; /**< dimensions */
+ int nel ; /**< number of image elements (pixels) */
+ int *subs ; /**< N-dimensional subscript */
+ int *dsubs ; /**< another subscript */
+ int *strides ; /**< strides to move in image data */
+ /*@}*/
+
+ vl_uint *perm ; /**< pixel ordering */
+ vl_uint *joins ; /**< sequence of join ops */
+ int njoins ; /**< number of join ops */
+
+ /** @name Regions */
+ /*@{*/
+ VlMserReg *r ; /**< basic regions */
+ VlMserExtrReg *er ; /**< extremal tree */
+ vl_uint *mer ; /**< maximally stable extremal regions */
+ int ner ; /**< number of extremal regions */
+ int nmer ; /**< number of maximally stable extr. reg. */
+ int rer ; /**< size of er buffer */
+ int rmer ; /**< size of mer buffer */
+ /*@}*/
+
+ /** @name Ellipsoids fitting */
+ /*@{*/
+ float *acc ; /**< moment accumulator. */
+ float *ell ; /**< ellipsoids list. */
+ int rell ; /**< size of ell buffer */
+ int nell ; /**< number of ellipsoids extracted */
+ int dof ; /**< number of dof of ellipsoids. */
+
+ /*@}*/
+
+ /** @name Configuration */
+ /*@{*/
+ vl_bool verbose ; /**< be verbose */
+ int delta ; /**< delta filter parameter */
+ double max_area ; /**< badness test parameter */
+ double min_area ; /**< badness test parameter */
+ double max_variation ; /**< badness test parameter */
+ double min_diversity ; /**< minimum diversity */
+ /*@}*/
+
+ VlMserStats stats ; /** run statistic */
+} ;
+
+/* ----------------------------------------------------------------- */
+/** @brief Get delta
+ ** @param f MSER filter.
+ ** @return value of @c delta.
+ **/
+VL_INLINE vl_mser_pix
+vl_mser_get_delta (VlMserFilt const *f)
+{
+ return f-> delta ;
+}
+
+/** @brief Set delta
+ ** @param f MSER filter.
+ ** @param x value of @c delta.
+ **/
+VL_INLINE void
+vl_mser_set_delta (VlMserFilt *f, vl_mser_pix x)
+{
+ f-> delta = x ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get minimum diversity
+ ** @param f MSER filter.
+ ** @return value of @c minimum diversity.
+ **/
+VL_INLINE double
+vl_mser_get_min_diversity (VlMserFilt const *f)
+{
+ return f-> min_diversity ;
+}
+
+/** @brief Set minimum diversity
+ ** @param f MSER filter.
+ ** @param x value of @c minimum diversity.
+ **/
+VL_INLINE void
+vl_mser_set_min_diversity (VlMserFilt *f, double x)
+{
+ f-> min_diversity = x ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get statistics
+ ** @param f MSER filter.
+ ** @return statistics.
+ **/
+VL_INLINE VlMserStats const*
+vl_mser_get_stats (VlMserFilt const *f)
+{
+ return & f-> stats ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get maximum region area
+ ** @param f MSER filter.
+ ** @return maximum region area.
+ **/
+VL_INLINE double
+vl_mser_get_max_area (VlMserFilt const *f)
+{
+ return f-> max_area ;
+}
+
+/** @brief Set maximum region area
+ ** @param f MSER filter.
+ ** @param x maximum region area.
+ **/
+VL_INLINE void
+vl_mser_set_max_area (VlMserFilt *f, double x)
+{
+ f-> max_area = x ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get minimum region area
+ ** @param f MSER filter.
+ ** @return minimum region area.
+ **/
+VL_INLINE double
+vl_mser_get_min_area (VlMserFilt const *f)
+{
+ return f-> min_area ;
+}
+
+/** @brief Set minimum region area
+ ** @param f MSER filter.
+ ** @param x minimum region area.
+ **/
+VL_INLINE void
+vl_mser_set_min_area (VlMserFilt *f, double x)
+{
+ f-> min_area = x ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get maximum region variation
+ ** @param f MSER filter.
+ ** @return maximum region variation.
+ **/
+VL_INLINE double
+vl_mser_get_max_variation (VlMserFilt const *f)
+{
+ return f-> max_variation ;
+}
+
+/** @brief Set maximum region variation
+ ** @param f MSER filter.
+ ** @param x maximum region variation.
+ **/
+VL_INLINE void
+vl_mser_set_max_variation (VlMserFilt *f, double x)
+{
+ f-> max_variation = x ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get maximally stable extremal regions
+ ** @param f MSER filter.
+ ** @return array of MSER pivots.
+ **/
+VL_INLINE vl_uint const *
+vl_mser_get_regions (VlMserFilt const* f)
+{
+ return f-> mer ;
+}
+
+/** @brief Get number of maximally stable extremal regions
+ ** @param f MSER filter.
+ ** @return number of MSERs.
+ **/
+VL_INLINE vl_uint
+vl_mser_get_regions_num (VlMserFilt const* f)
+{
+ return f-> nmer ;
+}
+
+/* ----------------------------------------------------------------- */
+/** @brief Get ellipsoids
+ ** @param f MSER filter.
+ ** @return ellipsoids.
+ **/
+VL_INLINE float const *
+vl_mser_get_ell (VlMserFilt const* f)
+{
+ return f-> ell ;
+}
+
+/** @brief Get number of degrees of freedom of ellipsoids
+ ** @param f MSER filter.
+ ** @return number of degrees of freedom.
+ **/
+VL_INLINE vl_uint
+vl_mser_get_ell_dof (VlMserFilt const* f)
+{
+ return f-> dof ;
+}
+
+/** @brief Get number of ellipsoids
+ ** @param f MSER filter.
+ ** @return number of ellipsoids
+ **/
+VL_INLINE vl_uint
+vl_mser_get_ell_num (VlMserFilt const* f)
+{
+ return f-> nell ;
+}
+
+/* VL_MSER */
+#endif
diff --git a/vl/sift.c b/vl/sift.c
new file mode 100644
index 0000000..b5c594a
--- /dev/null
+++ b/vl/sift.c
@@ -0,0 +1,2175 @@
+#if WIN32
+#pragma warning(disable:4244 4305)
+#endif
+
+/** @internal
+ ** @file sift.c
+ ** @brief Scale Invariant Feature Transform (SIFT) - Definition
+ ** @author Andrea Vedaldi
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+/**
+
+@file sift.h
+@brief Scale Invariant Feature Transform (SIFT)
+@author Andrea Vedaldi
+
+@par "Credits:" May people have contributed with suggestions and bug
+reports. Although the following list is certainly incomplete, we would
+like to thank: Brian Fulkerson, Wei Dong, Loic, Giuseppe, Liu, Erwin,
+P. Ivanov, and Q. S. Luo.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@section sift Scale Invariant Feature Transform
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+This library module implements a
+@ref sift-filter-usage "SIFT filter object",
+a reusable object to extract SIFT features from one or
+multiple images of the same size.
+
+- @ref sift-intro
+ - @ref sift-intro-detector
+ - @ref sift-intro-descriptor
+ - @ref sift-intro-extensions
+- @ref sift-usage
+- @ref sift-tech
+ - @ref sift-tech-ss
+ - @ref sift-tech-detector
+ - @ref sift-tech-detector-peak
+ - @ref sift-tech-detector-edge
+ - @ref sift-tech-detector-orientation
+ - @ref sift-tech-descriptor
+ - @ref sift-tech-descriptor-can
+ - @ref sift-tech-descriptor-image
+ - @ref sift-tech-descriptor-std
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@section sift-intro Overview
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+A SIFT feature is a selected image region (also called keypoint) with
+an associated descriptor. Keypoints are extracted by the <b>@ref
+sift-intro-detector "SIFT detector"</b> and their descriptors are
+computed by the <b>@ref sift-descriptor "SIFT descriptor"</b>. It is
+also common to use independently the SIFT detector (i.e. computing the
+keypoints without descriptors) or the SIFT descriptor (i.e. computing
+descriptors of custom keypoints).
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-intro-detector SIFT detector
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+@sa
+@ref sift-tech-ss "Scale space technical details",
+@ref sift-tech-detector "Detector technical details"
+
+A SIFT <em>keypoint</em> is a cricular image region with an
+orientation. It is described by a geometric <em>frame</em> of four
+parameters: the keypoint center coordinates @e x and @e y, its @e
+scale (the radius of the region), and its @e orientation (an angle
+expresed in radians). The SIFT detector uses as keypoints image
+structures which resemble &ldquo;blobs&rdquo;. By searching for blobs
+at multiple scales and positions, the SIFT detector is invariant (or,
+more accurately, covaraint) to translation, rotations, and rescaling
+of the image.
+
+The keypoint orientation is also determined from the local image
+appearance and is covariant to image rotations. Depending on the
+symmetry of the keypoint appearance, determining the orientation can
+be ambiguous. In this case, the SIFT detectors returns a list of up to
+four possible orientations, constructing up to four frames (differing
+only by their orientation) for each detected image blob.
+
+@image html sift-frame.png "SIFT keypoints are circular image regions with an orientation."
+
+There are several parameters that influence the detection of SIFT
+keypoints. First, searching keypoints at multiple scales is obtained
+by constructing a so-called &ldquo;Gausssian scale space&rdquo;. The
+scale space is just a collection of images obtained by progressively
+smoothing the input image, which is analogous to gradually reducing
+the image resolution. Conventionally, the smoothing level is called
+<em>scale</em> of the image. The construction of the scale space is
+influenced by the following parameters, set when creating the SIFT
+filter object by ::vl_sift_new():
+
+- <b>Number of octaves</b>. Increasing the scale by an octave means
+ doubling the size of the smoothing kernel, whose effect is roughly
+ equivalent to halving the image resolution. By default, the scale
+ space spans as many octaves as possible (i.e. roughy <code>
+ log2(min(width,height)</code>), which has the effect of searching
+ keypoints of all possible sizes.
+- <b>First octave index</b>. By convention, the octave of index 0
+ starts with the image full resolution. Specifying an index greater
+ than 0 starts the scale space at a lower resolution (e.g. 1 halves
+ the resolution). Similarly, specifying a negative index starts the
+ scale space at an higher resoltuon image, and can be useful to
+ extract very small features (since this is obtained by interpolating
+ the input image, it does not make much sense to go past -1).
+- <b>Number of levels per octave</b>. Each octave is sampled at this
+ given numer of intermediate scales (by default 3). Increasing this
+ number might in principle return more refined keypoints, but in
+ practice can make their selection unstable due to noise (see [1]).
+
+Keypoints are further refined by eliminating those that are likely to
+be unstable, either because they are selected nearby an image edge,
+rather than an image blob, or are found on image structures with low
+contrast. Filtering is controlled by the follow:
+
+- <b>Peak threshold.</b> This is the minimum amount of contrast to
+ accept a keypoint. It is set by configuring the SIFT filter object
+ by ::vl_sift_set_peak_thresh().
+- <b>Edge threhsold.</b> This is the edge rejection threhsold. It is
+ set by configuring the SIFT filter object by
+ ::vl_sift_set_edge_thresh().
+
+<table>
+ <caption>Summary of the parameters influencing the SIFT detector.</caption>
+ <tr style="font-weight:bold;">
+ <td>Parameter</td>
+ <td>See also</td>
+ <td>Controlled by</td>
+ <td>Comment</td>
+ </tr>
+ <tr>
+ <td>number of octaves</td>
+ <td> @ref sift-intro-detector </td>
+ <td>::vl_sift_new</td>
+ <td></td>
+ </tr>
+ <tr>
+ <td>first octave index</td>
+ <td> @ref sift-intro-detector </td>
+ <td>::vl_sift_new</td>
+ <td>set to -1 to extract very small features</td>
+ </tr>
+ <tr>
+ <td>number of scale levels per octave</td>
+ <td> @ref sift-intro-detector </td>
+ <td>::vl_sift_new</td>
+ <td>can affect the number of extracted keypoints</td>
+ </tr>
+ <tr>
+ <td>edge threhsold</td>
+ <td> @ref sift-intro-detector </td>
+ <td>::vl_sift_set_edge_thresh</td>
+ <td>decrease to eliminate more keypoints</td>
+ </tr>
+ <tr>
+ <td>peak threhsold</td>
+ <td> @ref sift-intro-detector </td>
+ <td>::vl_sift_set_peak_thresh</td>
+ <td>increase to eliminate more keypoints</td>
+ </tr>
+</table>
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-intro-descriptor SIFT Descriptor
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+@sa @ref sift-tech-descriptor "Descriptor technical details"
+
+A SIFT descriptor is a 3-D spatial histogram of the image gradients in
+characterizing the appearance of a keypoint. The gradient at each
+pixel is regarded as a sample of a three-dimensional elementary
+feature vector, formed by the pixel location and the gradient
+orientation. Samples are weighed by the gradient norm and accumulated
+in a 3-D histogram @em h, which (up to normalization and clamping)
+forms the SIFT descriptor of the region. An additional Gaussian
+weighting function is applied to give less importance to gradients
+farther away from the keypoint center. Orientations are quantized into
+eight bins and the spatial coordinates into four each, as follows:
+
+@image html sift-descr-easy.png "The SIFT descriptor is a spatial histogram of the image gradient."
+
+SIFT descriptors are computed by either calling
+::vl_sift_calc_keypoint_descriptor or
+::vl_sift_cal_keypoint_descriptor_raw. They accept as input a keypoint
+frame, which specifies the descriptor center, its size, and its
+orientation on the image plane. The following parameters influence the
+descriptor calculation:
+
+- <b>magnification factor</b>. The desciriptor size is determined by
+multiplying the keypoint scale by this factor. It is set by
+::vl_sift_set_magnif.
+- <b>Gaussian window size</b>. The desriptor support is determined by
+a Gaussian window, which discounts gradient contributions farther away
+from the descriptor center. The standard deviation of this window is
+set by ::vl_sift_set_window_size and expressed in unit of bins.
+
+VLFeat SIFT descriptor uses the following convention. The @em y axis
+points downards and angles are measured clockwise (to be consistent
+with the standard image convention). The 3-D histogram (consisting of
+@f$ 8 \times 4 \times 4 = 128 @f$ bins) is stacked as a single
+128-dimensional vector, where the fastest varying dimension is the
+orientation and the slowest the @em y spatial coordinate. This is
+illustrated by the following figure.
+
+@image html sift-conv-vlfeat.png "VLFeat conventions"
+
+@note Keypoints (frames) D. Lowe's SIFT implementation convention is
+slightly different: The @em y axis points upwards and the angles are
+measured counter-clockwise.
+
+@image html sift-conv.png "D. Lowes' SIFT implementation conventions"
+
+<table>
+ <caption>Summary of the parameters influencing the SIFT descriptor.</caption>
+ <tr style="font-weight:bold;">
+ <td>Parameter</td>
+ <td>See also</td>
+ <td>Controlled by</td>
+ <td>Comment</td>
+ </tr>
+ <tr>
+ <td>magnification factor</td>
+ <td> @ref sift-intro-descriptor </td>
+ <td>::vl_sift_set_magnif</td>
+ <td>increase this value to enlarge the image region described</td>
+ </tr>
+ <tr>
+ <td>Gaussian window size</td>
+ <td> @ref sift-intro-descriptor </td>
+ <td>::vl_sift_set_window_size</td>
+ <td>smaller values let the center of the descriptor count more</td>
+ </tr>
+</table>
+
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@section sift-intro-extensions Extensions
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+<b>Eliminating low-contrast descriptors.</b> Near-uniform patches do
+not yield stable keypoints or descriptors. ::vl_sift_set_norm_thresh()
+can be used to set a threshold on the average norm of the local
+gradient to zero-out descriptors that correspond to very low contrast
+regions. By default, the threshold is equal to zero, which means that
+no descriptor is zeroed. Normally this option is useful only with
+custom keypoints, as detected keypoints are implicitly eseleted at
+high contrast image regions.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@section sift-usage Using the SIFT filter object
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+The code provided in this module can be used in different ways. You
+can instantiate and use a <b>SIFT filter</b> to extract both SIFT
+keypoints and descriptors from one or multiple images. Alternatively,
+you can use one of the low level functions to run only a part of the
+SIFT algorithm (for instance, to compute the SIFT descriptors of
+custom keypoints).
+
+To use a <b>SIFT filter</b> object:
+
+- Initialize a SIFT filter object with ::vl_sift_new(). The filter can
+ be reused for multiple images of the same size (e.g. for an entiere
+ video sequence).
+- For each octave in the scale space:
+ - Compute the next octave of the DOG scale space using either
+ ::vl_sift_process_first_octave() or ::vl_sift_process_next_octave()
+ (stop processing if ::VL_ERR_EOF is returned).
+ - Run the SIFT detector with ::vl_sift_detect() to get the keypoints.
+ - For each keypoint:
+ - Use ::vl_sift_calc_keypoint_orientations() to get the keypoint orientation(s).
+ - For each orientation:
+ - Use ::vl_sift_calc_keypoint_descriptor() to get the keypoint descriptor.
+- Delete the SIFT filter by ::vl_sift_delete().
+
+To compute SIFT descriptors of custom keypoints, use
+::vl_sift_calc_raw_descriptor().
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@section sift-tech Technical details
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-tech-ss Scale space
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+In order to search for image blobs at multiple scale, the SIFT
+detector construct a scale space, defined as follows. Let
+@f$I_0(\mathbf{x})@f$ denote an idealized <em>infinite resolution</em>
+image. Consider the <em>Gaussian kernel</em>
+
+@f[
+ g_{\sigma}(\mathbf{x})
+ =
+ \frac{1}{2\pi\sigma^2}
+ \exp
+ \left(
+ -\frac{1}{2}
+ \frac{\mathbf{x}^\top\mathbf{x}}{\sigma^2}
+ \right)
+@f]
+
+The <b>Gaussian scale space</b> is the collection of smoothed images
+
+@f[
+ I_\sigma = g_\sigma * I, \quad \sigma \geq 0.
+@f]
+
+The image at infinite resolution @f$ I_0 @f$ is useful conceptually,
+but is not available to us; intstead, the input image @f$ I_{\sigma_n}
+@f$ is assumed to be pre-smoothed at a nomimal level @f$ \sigma_n =
+0.5 @f$ to account for the finite resolution of the pixels. Thus in
+practice the scale space is computed by
+
+@f[
+I_\sigma = g_{\sqrt{\sigma^2 - \sigma_n^2}} * I_{\sigma_n},
+\quad \sigma \geq \sigma_n.
+@f]
+
+Scales are sampled at logarithmic steps given by
+
+@f[
+\sigma = \sigma_0 2^{o+s/S},
+\quad s = 0,\dots,S-1,
+\quad o = o_{\min}, \dots, o_{\min}+O-1,
+@f]
+
+being @f$ \sigma_0 = 1.6 @f$ is the <em>base scale</em>, @f$ o_{\min}
+@f$ is the <em>first octave index</em>, @em O the <em>number of
+octaves</em> and @em S the <em>number of scales per octave</em>.
+
+Blobs are detected as local extrema of the <b>Difference of
+Gaussians</b> (DoG) scale space, obtained by subtracting successive
+scales of the Gaussian scale space:
+
+@f[
+\mathrm{DoG}_{\sigma(o,s)} = I_{\sigma(o,s+1)} - I_{\sigma(o,s)}
+@f]
+
+At each next octave, the resoltuion of the images is halved to save
+computations. The images composing the Gaussian and DoG scale space
+can then be arranged as in the following figure:
+
+@image html sift-ss.png "GSS and DoG scale space structures."
+
+The black vertical segments represent images of the Gaussian Scale
+Space (GSS), arranged by increasing scale @f$\sigma@f$. Notice that
+the sacle level index @e s varies in a slightly redundant set
+
+@f[
+s = -1, \dots, S+2
+@f]
+
+This simplifies gluging togheter different octaves and extracting DoG
+maxima (required by the SIFT detector).
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-tech-detector Detector
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+The SIFT frames (keypoints) are extracted based on local extrema
+(peaks) of the DoG scale space. Numerically, local extrema are
+elements whose @f$ 3 \times 3 \times 3 @f$ neighbors (in space and
+scale) have all smaller (or larger) value. Once extracted, local
+extrema are quadratically interpolated (this is very important
+especially at the lower resolution scales in order to have accurate
+keypoint localization at the full resolution). Finally, they are
+filtered to eliminate low-contranst responses or responses close to
+edges and the orientation(s) are assigned, as explained next.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsubsection sift-tech-detector-peak Eliminating low contrast responses
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+Peaks which are too short may have been generated by noise and are
+discarded. This is done by comparing the absolute value of the DoG
+scale space at the peak with the <b>peak threshold</b> @f$t_p@f$ and
+discarding the peak its value is below the threshold.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsubsection sift-tech-detector-edge Eliminating edge responses
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+Peaks which are too flat are often generated by edges and do not yield
+stable features. These peaks are detected and removed as follows.
+Given a peak @f$x,y,\sigma@f$, the algorithm evaluates the @em x,@em y
+Hessian of of the DoG scale space at the scale @f$\sigma@f$. Then the
+following score (similar to the Harris function) is computed:
+
+@f[
+\frac{(\mathrm{tr}\,D(x,y,\sigma))^2}{\det D(x,y,\sigma)},
+\quad
+D =
+\left[
+\begin{array}{cc}
+\frac{\partial^2 \mathrm{DoG}}{\partial x^2} &
+\frac{\partial^2 \mathrm{DoG}}{\partial x\partial y} \\
+\frac{\partial^2 \mathrm{DoG}}{\partial x\partial y} &
+\frac{\partial^2 \mathrm{DoG}}{\partial y^2}
+\end{array}
+\right].
+ @f]
+
+This score has a minimum (equal to 4) when both eigenvalues of the
+Jacobian are equal (curved peak) and increases as one of the
+eigenvalues grows and the other stays small. Peaks are retained if the
+score is below the quantity @f$(t_e+1)(t_e+1)/t_e@f$, where @f$t_e@f$
+is the <b>edge threshold</b>. Notice that this quantity has a minimum
+equal to 4 when @f$t_e=1@f$ and grows thereafter. Therefore the range
+of the edge threshold is @f$[1,\infty)@f$.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-tech-detector-orientation Orientation assignment
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+A peak in the DoG scale space fixes 2 parameters of the keypoint: the
+position and scale. It remains to choose an orientation. In order to
+do this, SIFT computes an histogram of the gradient orientations in a
+Gaussian window with a standard deviation which is 1.5 times bigger
+than the scale @f$\sigma@f$ of the keypoint.
+
+@image html sift-orient.png
+
+This histogram is then smoothed and the maximum is selected. In
+addition to the biggest mode, up to other three modes whose amplitude
+is within the 80% of the biggest mode are retained and returned as
+additional orientations.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsection sift-tech-descriptor Descriptor
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+A SIFT descriptor of a local region (keypoint) is a 3-D spatial
+histogram of the image gradients. The gradient at each pixel is
+regarded as a sample of a three-dimensional elementary feature vector,
+formed by the pixel location and the gradient orientation. Samples are
+weighed by the gradient norm and accumulated in a 3-D histogram @em h,
+which (up to normalization and clamping) forms the SIFT descriptor of
+the region. An additional Gaussian weighting function is applied to
+give less importance to gradients farther away from the keypoint
+center.
+
+<!-- @image html sift-bins.png -->
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsubsection sift-tech-descriptor-can Construction in the canonical frame
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+Denote the gradient vector field computed at the scale @f$ \sigma @f$ by
+@f[
+ J(x,y) = \nalba I_\sigma(x,y)
+ =
+ \left[\begin{array}{cc}
+ \frac{\partial I_\sigma}{\partial x} &
+ \frac{\partial I_\sigma}{\partial y} &
+ \end{array}\right]
+@f]
+
+The descriptor is a 3-D spatial histogram capturing the distribution
+of @f$ J(x,y) @f$. It is convenient to describe its construction in
+the <em>canonical frame</em>. In this frame, the image and descriptor
+axes coincide and each spatial bin has side 1. The histogram has @f$
+N_\theta \times N_x \times N_y @f$ bins (usually @f$ 8 \times 4 \times
+4 @f$), as in the following figure:
+
+@image html sift-can.png Canonical SIFT descriptor and spatial binning functions
+
+Bins are indexed by a triplet of indexes <em>t, i, j</em> and their
+centers are given by
+
+@f{eqnarray*}
+ \theta_t &=& \frac{2\pi}{N_\theta} t, \quad t = 0,\dots,N_{\theta}-1, \\
+ x_i &=& i - \frac{N_x -1}{2}, \quad i = 0,\dots,N_x-1, \\
+ y_j &=& j - \frac{N_x -1}{2}, \quad j = 0,\dots,N_y-1. \\
+@f}
+
+The histogram is computed by using trilinear interpolation, i.e. by
+weighing contributions by the <em>binning functions</em>
+
+@f{eqnarray*}
+ \displaystyle
+ w(z) &=& \mathrm{max}(0, 1 - |z|),
+ \\
+ \displaystyle
+ w_\mathrm{ang}(z) &=& \sum_{k=-\inf}^{+\inf}
+ w\left(
+ \frac{N_\theta}{2\pi} z + N_\theta k
+ \right).
+@f}
+
+The gradient vector field is transformed in a three-dimensional
+density map of weighed contribtions
+
+@f[
+ f(\theta, x, y) = |J(x,y)| \delta(\theta - \angle J(x,y))
+@f]
+
+The historam is localized in the keypoint support by a Gaussian window
+of standard deviation @f$ \sigma_{\mathrm{win}} @f$. The histogram is
+then given by
+
+@f{eqnarray*}
+ h(t,i,j) &=& \int
+ g_{\sigma_\mathrm{win}}(x,y)
+ w_\mathrm{ang}(\theta - \theta_t) w(x-x_i) w(y-y_j)
+ f(\theta,x,y)
+ d\theta\,dx\,dy
+\\
+&=& \int
+ g_{\sigma_\mathrm{win}}(x,y)
+ w_\mathrm{ang}(\angle J(x,y) - \theta_t) w(x-x_i) w(y-y_j)
+ |J(x,y)|\,dx\,dy
+@f}
+
+In post processing, the histogram is @f$ l^2 @f$ normalized, then
+clamped at 0.2, and @f$ l^2 @f$ normalized again.
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsubsection sift-tech-descriptor-image Calculation in the image frame
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+Invariance to similarity transformation is attained by atttaching
+descriptors to SIFT keypoints (or other similarity-covariant frames).
+Then projecting the image in the canonical descriptor frames has the
+effect of undoing the image deformation.
+
+In practice, however, it is convenient to compute the descriptor
+directly in the image frame. To do this, denote with a hat quantities
+relative to the canonical frame and without a hat quantities relative
+to the image frame (so for instance @f$ \hat x @f$ is the @e
+x-coordinate in the canonical frame and @f$ x @f$ the x-coordinate in
+the image frame). Assume that canonical and image frame are
+related by an affinity:
+
+@f[
+ \mathbf{x} = A \hat\mathbf{x} + T,
+ \qquad
+ \mathbf{x} =
+ \left[\begin{array}{cc}
+ x \\
+ y
+ \end{arraty}\right],
+ \quad
+ \hat\mathbf{x} =
+ \left[\begin{array}{cc}
+ \hat x \\
+ \hat y
+ \end{arraty}\right].
+@f]
+
+@image html sift-image-frame.png
+
+Then all quantites can be computed in the image frame directly. For instance,
+the image at infinite resolution in the two frames are related by
+
+@f[
+ \hat I_0(\hat\mathbf{x}) = I_0(\mathbf{x}), \qquad \mathbf{x} = A \hat\mathbf{x} + T.
+@f]
+
+The canonized image at scale @f$ \hat \sigma @f$ is in relation with the scaled image
+
+@f[
+ \hat I_{\hat{\sigma}}(\hat\mathbf{x}) = I_{A\hat{\sigma}}(\mathbf{x}),
+ \qquad \mathbf{x} = A \hat\mathbf{x} + T
+@f]
+
+where by generalizing the previous definitions we have
+
+@f[
+ I_{A\hat \sigma}(\mathbf{x}) = (g_{A\hat\sigma} * I_0)(\mathbf{x}),
+\quad
+ g_{A\hat\sigma}(\mathbf{x})
+ =
+ \frac{1}{2\pi|A|\hat \sigma^2}
+ \exp
+ \left(
+ -\frac{1}{2}
+ \frac{\mathbf{x}^\top A^{-\top}A^{-1}\mathbf{x}}{\hat \sigma^2}
+ \right)
+@f]
+
+Deriving shows that the gradient fileds are in relation
+
+@f[
+ \hat J(\hat \mathbf{x}) = J(\mathbf{x}) A,
+ \quad J(\mathbf{x}) = (\nabla I_{A\hat\sigma})(\mathbf{x}),
+ \qquad \mathbf{x} = A \hat\mathbf{x} + T.
+@f]
+
+So we can compute the descriptor either in the image or canonical frame as:
+
+@f{eqnarray*}
+ h(t,i,j)
+ &=&
+ \int
+ g_{\hat \sigma_\mathrm{win}}(\hat \mathbf{x})\,
+ w_\mathrm{ang}(\angle \hat J(\hat\mathbf{x}) - \theta_t)\,
+ w_{ij}(\hat\mathbf{x})\,
+ |\hat J(\hat \mathbf{x})|\,
+ d\hat \mathbf{x}
+ \\
+ &=& \int
+ g_{A \hat \sigma_\mathrm{win}}(\mathbf{x} - T)\,
+ w_\mathrm{ang}(\angle J(\mathbf{x})A - \theta_t)\,
+ w_{ij}(A^{-1}(\mathbf{x} - T))\,
+ |J(\mathbf{x})A|\,
+ d\mathbf{x}.
+@f}
+
+where we defined the product of the two spatial binning functions
+
+@f[
+ w_{ij}(\hat\mathbf{x}) = w(\hat x - \hat x_i) w(\hat y - \hat y_j)
+@f]
+
+
+In the actual implementation, this integral is computed by visiting a
+rectangular area of the image that fully contains the keypoint grid
+(along with half a bin border to fully include the bin windowing
+function). Since the descriptor can be rotated, this area is a
+rectangle of sides @f$m/2\sqrt{2} (N_x+1,N_y+1)@f$ (see also the
+illustration).
+
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+@subsubsection sift-tech-descriptor-std Standard SIFT descriptor
+<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
+
+For a SIFT-detected keypoint of center @f$ T @f$, scale @f$ \sigma @f$
+and orientation @f$ \theta @f$, the affine transformation @f$ (A,T)
+@f$ reduces to the similatity tranformation
+
+@f[
+ \mathbf{x} = m \sigma R(\theta) \hat \mathbf{x} + T
+@f]
+
+where @f$ R(\theta) @f$ is a counter-clockwise rotation of @f$ \theta
+@f$ radians, and @e m is the <b>descriptor magnification factor</b>
+which expresses how much larger the descriptor window is compared to
+the scale of the keypoint (a common value is @e m = 3). Moreover, the
+standard SIFT descriptor computes the image gradient at the scale of
+the keypoints, which in the canonical frame is equivalent to a
+smoothing of @f$ \hat \sigma = 1/m @f$. Finally, the Gaussian window
+is set to have standard deviation @f$ \hat \sigma_\mathrm{win} = 2
+@f$. This yields the formula
+
+@f{eqnarray*}
+ h(t,i,j)
+ &=&
+ m \sigma \int
+ g_{\sigma_\mathrm{win}}(\mathbf{x} - T)\,
+ w_\mathrm{ang}(\angle J(\mathbf{x}) - \theta - \theta_t)\,
+ w_{ij}\left(\frac{R(\theta)^\top \mathbf{x} - T}{m\sigma}\right)\,
+ |J(\mathbf{x})|\,
+ d\mathbf{x},
+\\
+\sigma_{\mathrm{win}} &=& m\sigma\hat \sigma_{\mathrm{win}},
+\\
+ J(\mathbf{x})
+ &=& \nabla (g_{m \sigma \hat \sigma} * I)(\mathbf{x})
+ = \nabla (g_{\sigma} * I)(\mathbf{x})
+ = \nabla I_{\sigma} (\mathbf{x}).
+@f}
+
+
+
+**/
+
+#include "sift.h"
+#include "imop.h"
+#include "mathop.h"
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdio.h>
+
+/** @internal @brief Use bilinear interpolation to compute orientations */
+#define VL_SIFT_BILINEAR_ORIENTATIONS 1
+
+#define EXPN_SZ 256 /**< ::fast_expn table size @internal */
+#define EXPN_MAX 25.0 /**< ::fast_expn table max @internal */
+double expn_tab [EXPN_SZ+1] ; /**< ::fast_expn table @internal */
+
+#define NBO 8
+#define NBP 4
+
+#define log2(x) (log(x)/VL_LOG_OF_2)
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Fast @f$exp(-x)@f$ approximation
+ **
+ ** @param x argument.
+ **
+ ** The argument must be in the range [0, ::EXPN_MAX] .
+ **
+ ** @return approximation of @f$exp(-x)@f$.
+ **/
+
+VL_INLINE double
+fast_expn (double x)
+{
+ double a,b,r ;
+ int i ;
+ /*assert(0 <= x && x <= EXPN_MAX) ;*/
+
+ if (x > EXPN_MAX) return 0.0 ;
+
+ x *= EXPN_SZ / EXPN_MAX ;
+ i = vl_floor_d (x) ;
+ r = x - i ;
+ a = expn_tab [i ] ;
+ b = expn_tab [i + 1] ;
+ return a + r * (b - a) ;
+}
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Initialize tables for ::fast_expn
+ **/
+
+VL_INLINE void
+fast_expn_init ()
+{
+ int k ;
+ for(k = 0 ; k < EXPN_SZ + 1 ; ++ k) {
+ expn_tab [k] = exp (- (double) k * (EXPN_MAX / EXPN_SZ)) ;
+ }
+}
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Copy image, upsample rows and take transpose
+ **
+ ** @param dst output image buffer.
+ ** @param src input image buffer.
+ ** @param width input image width.
+ ** @param height input image height.
+ **
+ ** The output image has dimensions @a height by 2 @a width (so the
+ ** destination buffer must be at least as big as two times the
+ ** input buffer).
+ **
+ ** Upsampling is performed by linear interpolation.
+ **/
+
+static void
+copy_and_upsample_rows
+(vl_sift_pix *dst,
+ vl_sift_pix const *src, int width, int height)
+{
+ int x, y ;
+ vl_sift_pix a, b ;
+
+ for(y = 0 ; y < height ; ++y) {
+ b = a = *src++ ;
+ for(x = 0 ; x < width - 1 ; ++x) {
+ b = *src++ ;
+ *dst = a ; dst += height ;
+ *dst = 0.5 * (a + b) ; dst += height ;
+ a = b ;
+ }
+ *dst = b ; dst += height ;
+ *dst = b ; dst += height ;
+ dst += 1 - width * 2 * height ;
+ }
+}
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Copy and downsample an image
+ **
+ ** @param dst output imgae buffer.
+ ** @param src input image buffer.
+ ** @param width input image width.
+ ** @param height input image height.
+ ** @param d octaves (non negative).
+ **
+ ** The function downsamples the image @a d times, reducing it to @c
+ ** 1/2^d of its original size. The parameters @a width and @a height
+ ** are the size of the input image. The destination image @a dst is
+ ** assumed to be <code>floor(width/2^d)</code> pixels wide and
+ ** <code>floor(height/2^d)</code> pixels high.
+ **/
+
+static void
+copy_and_downsample
+(vl_sift_pix *dst,
+ vl_sift_pix const *src,
+ int width, int height, int d)
+{
+ int x, y ;
+
+ d = 1 << d ; /* d = 2^d */
+ for(y = 0 ; y < height ; y+=d) {
+ vl_sift_pix const * srcrowp = src + y * width ;
+ for(x = 0 ; x < width - (d-1) ; x+=d) {
+ *dst++ = *srcrowp ;
+ srcrowp += d ;
+ }
+ }
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Create a new SIFT filter
+ **
+ ** @param width image width.
+ ** @param height image height.
+ ** @param noctaves number of octaves.
+ ** @param nlevels number of levels per octave.
+ ** @param o_min first octave index.
+ **
+ ** The function allocates and returns a new SIFT filter for the
+ ** specified image and scale space geometry.
+ **
+ ** Setting @a O to a negative value sets the number of octaves to the
+ ** maximum possible value depending on the size of the image.
+ **
+ ** @return the new SIFT filter.
+ ** @sa ::vl_sift_delete().
+ **/
+
+VL_EXPORT
+VlSiftFilt *
+vl_sift_new (int width, int height,
+ int noctaves, int nlevels,
+ int o_min)
+{
+ VlSiftFilt *f = vl_malloc (sizeof(VlSiftFilt)) ;
+
+ int w = VL_SHIFT_LEFT (width, -o_min) ;
+ int h = VL_SHIFT_LEFT (height, -o_min) ;
+ int nel = w * h ;
+
+ /* negative value O => calculate max. value */
+ if (noctaves < 0) {
+ noctaves = VL_MAX (floor (log2 (VL_MIN(width, height))) - o_min - 3, 1) ;
+ }
+
+ f-> width = width ;
+ f-> height = height ;
+ f-> O = noctaves ;
+ f-> S = nlevels ;
+ f-> o_min = o_min ;
+ f-> s_min = -1 ;
+ f-> s_max = nlevels + 1 ;
+ f-> o_cur = o_min ;
+
+ f-> temp = vl_malloc (sizeof(vl_sift_pix) * nel ) ;
+ f-> octave = vl_malloc (sizeof(vl_sift_pix) * nel
+ * (f->s_max - f->s_min + 1) ) ;
+ f-> dog = vl_malloc (sizeof(vl_sift_pix) * nel
+ * (f->s_max - f->s_min ) ) ;
+ f-> grad = vl_malloc (sizeof(vl_sift_pix) * nel * 2
+ * (f->s_max - f->s_min ) ) ;
+
+ f-> sigman = 0.5 ;
+ f-> sigma0 = 1.6 * pow (2.0, 1.0 / nlevels) ;
+ f-> sigmak = pow (2.0, 1.0 / nlevels) ;
+ f-> dsigma0 = f->sigma0 * sqrt (1.0 - 1.0 / (f->sigmak*f->sigmak)) ;
+
+ /*
+ VL_PRINTF ("sigman = %g\n", f-> sigman) ;
+ VL_PRINTF ("sigma0 = %g\n", f-> sigma0) ;
+ */
+
+ f-> octave_width = 0 ;
+ f-> octave_height = 0 ;
+
+ f-> keys = 0 ;
+ f-> nkeys = 0 ;
+ f-> keys_res = 0 ;
+
+ f-> peak_thresh = 0.0 ;
+ f-> edge_thresh = 10.0 ;
+ f-> norm_thresh = 0.0 ;
+ f-> magnif = 3.0 ;
+ f-> windowSize = NBP / 2 ;
+
+ f-> grad_o = o_min - 1 ;
+
+ /* initialize fast_expn stuff */
+ fast_expn_init () ;
+
+ return f ;
+}
+
+/** -------------------------------------------------------------------
+ ** @brief Delete SIFT filter
+ **
+ ** @param f SIFT filter to delete.
+ **
+ ** The function frees the resources allocated by ::vl_sift_new().
+ **/
+
+VL_EXPORT
+void
+vl_sift_delete (VlSiftFilt* f)
+{
+ if(f) {
+ if(f-> keys ) vl_free (f-> keys ) ;
+ if(f-> grad ) vl_free (f-> grad ) ;
+ if(f-> dog ) vl_free (f-> dog ) ;
+ if(f-> octave ) vl_free (f-> octave ) ;
+ if(f-> temp ) vl_free (f-> temp ) ;
+ vl_free (f) ;
+ }
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Start processing a new image
+ **
+ ** @param f SIFT filter.
+ ** @param im image data.
+ **
+ ** The function starts processing a new image by computing its
+ ** Gaussian scale space at the lower octave. It also empties the
+ ** internal keypoint buffer.
+ **
+ ** @return error code. The function returns ::VL_ERR_EOF if there are
+ ** no more octaves to process.
+ **
+ ** @sa ::vl_sift_process_next_octave().
+ **/
+
+VL_EXPORT
+int
+vl_sift_process_first_octave (VlSiftFilt *f, vl_sift_pix const *im)
+{
+ int o, s, h, w ;
+ double sa, sb ;
+ vl_sift_pix *octave ;
+
+ /* shortcuts */
+ vl_sift_pix *temp = f-> temp ;
+ int width = f-> width ;
+ int height = f-> height ;
+ int o_min = f-> o_min ;
+ int s_min = f-> s_min ;
+ int s_max = f-> s_max ;
+ double sigma0 = f-> sigma0 ;
+ double sigmak = f-> sigmak ;
+ double sigman = f-> sigman ;
+ double dsigma0 = f-> dsigma0 ;
+
+ /* restart from the first */
+ f->o_cur = o_min ;
+ f->nkeys = 0 ;
+ w = f-> octave_width = VL_SHIFT_LEFT(f->width, - f->o_cur) ;
+ h = f-> octave_height = VL_SHIFT_LEFT(f->height, - f->o_cur) ;
+
+ /* is there at least one octave? */
+ if (f->O == 0)
+ return VL_ERR_EOF ;
+
+ /* ------------------------------------------------------------------
+ * Compute the first sublevel of the first octave
+ * --------------------------------------------------------------- */
+
+ /*
+ * If the first octave has negative index, we upscale the image; if
+ * the first octave has positive index, we downscale the image; if
+ * the first octave has index zero, we just copy the image.
+ */
+
+ octave = vl_sift_get_octave (f, s_min) ;
+
+ if (o_min < 0) {
+ /* double once */
+ copy_and_upsample_rows (temp, im, width, height) ;
+ copy_and_upsample_rows (octave, temp, height, 2 * width ) ;
+
+ /* double more */
+ for(o = -1 ; o > o_min ; --o) {
+ copy_and_upsample_rows (temp, octave,
+ width << -o, height << -o ) ;
+ copy_and_upsample_rows (octave, temp,
+ width << -o, 2 * (height << -o)) ;
+ }
+ }
+ else if (o_min > 0) {
+ /* downsample */
+ copy_and_downsample (octave, im, width, height, o_min) ;
+ }
+ else {
+ /* direct copy */
+ memcpy(octave, im, sizeof(vl_sift_pix) * width * height) ;
+ }
+
+ /*
+ * Here we adjust the smoothing of the first level of the octave.
+ * The input image is assumed to have nominal smoothing equal to
+ * f->simgan.
+ */
+
+ sa = sigma0 * pow (sigmak, s_min) ;
+ sb = sigman * pow (2.0, - o_min) ;
+
+ if (sa > sb) {
+ double sd = sqrt (sa*sa - sb*sb) ;
+ vl_imsmooth_f (octave, temp, octave, w, h, sd) ;
+ }
+
+ /* -----------------------------------------------------------------
+ * Compute the first octave
+ * -------------------------------------------------------------- */
+
+ for(s = s_min + 1 ; s <= s_max ; ++s) {
+ double sd = dsigma0 * pow (sigmak, s) ;
+ vl_imsmooth_f (vl_sift_get_octave(f, s ), temp,
+ vl_sift_get_octave(f, s - 1), w, h, sd) ;
+ }
+
+ return VL_ERR_OK ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Process next octave
+ **
+ ** @param f SIFT filter.
+ **
+ ** The function computes the next octave of the Gaussian scale space.
+ ** Notice that this clears the record of any feature detected in the
+ ** previous octave.
+ **
+ ** @return error code. The function returns the error
+ ** ::VL_ERR_EOF when there are no more octaves to process.
+ **
+ ** @sa ::vl_sift_process_first_octave().
+ **/
+
+VL_EXPORT
+int
+vl_sift_process_next_octave (VlSiftFilt *f)
+{
+
+ int s, h, w, s_best ;
+ double sa, sb ;
+ vl_sift_pix *octave, *pt ;
+
+ /* shortcuts */
+ vl_sift_pix *temp = f-> temp ;
+ int O = f-> O ;
+ int S = f-> S ;
+ int o_min = f-> o_min ;
+ int s_min = f-> s_min ;
+ int s_max = f-> s_max ;
+ double sigma0 = f-> sigma0 ;
+ double sigmak = f-> sigmak ;
+ double dsigma0 = f-> dsigma0 ;
+
+ /* is there another octave ? */
+ if (f->o_cur == o_min + O - 1)
+ return VL_ERR_EOF ;
+
+ /* retrieve base */
+ s_best = VL_MIN(s_min + S, s_max) ;
+ w = vl_sift_get_octave_width (f) ;
+ h = vl_sift_get_octave_height (f) ;
+ pt = vl_sift_get_octave (f, s_best) ;
+ octave = vl_sift_get_octave (f, s_min) ;
+
+ /* next octave */
+ copy_and_downsample (octave, pt, w, h, 1) ;
+
+ f-> o_cur += 1 ;
+ f-> nkeys = 0 ;
+ w = f-> octave_width = VL_SHIFT_LEFT(f->width, - f->o_cur) ;
+ h = f-> octave_height = VL_SHIFT_LEFT(f->height, - f->o_cur) ;
+
+ sa = sigma0 * powf (sigmak, s_min ) ;
+ sb = sigma0 * powf (sigmak, s_best - S) ;
+
+ if (sa > sb) {
+ double sd = sqrt (sa*sa - sb*sb) ;
+ vl_imsmooth_f (octave, temp, octave, w, h, sd) ;
+ }
+
+ /* ------------------------------------------------------------------
+ * Fill octave
+ * --------------------------------------------------------------- */
+
+ for(s = s_min + 1 ; s <= s_max ; ++s) {
+ double sd = dsigma0 * pow (sigmak, s) ;
+ vl_imsmooth_f (vl_sift_get_octave(f, s ), temp,
+ vl_sift_get_octave(f, s - 1), w, h, sd) ;
+ }
+
+ return VL_ERR_OK ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Detect keypoints
+ **
+ ** The function detect keypoints in the current octave filling the
+ ** internal keypoint buffer. Keypoints can be retrieved by
+ ** ::vl_sift_get_keypoints().
+ **
+ ** @param f SIFT filter.
+ **/
+
+VL_EXPORT
+void
+vl_sift_detect (VlSiftFilt * f)
+{
+ vl_sift_pix* dog = f-> dog ;
+ int s_min = f-> s_min ;
+ int s_max = f-> s_max ;
+ int w = f-> octave_width ;
+ int h = f-> octave_height ;
+ double te = f-> edge_thresh ;
+ double tp = f-> peak_thresh ;
+
+ int const xo = 1 ; /* x-stride */
+ int const yo = w ; /* y-stride */
+ int const so = w * h ; /* s-stride */
+
+ double xper = pow (2.0, f->o_cur) ;
+
+ int x, y, s, i, ii, jj ;
+ vl_sift_pix *pt, v ;
+ VlSiftKeypoint *k ;
+
+ /* clear current list */
+ f-> nkeys = 0 ;
+
+ /* compute difference of gaussian (DoG) */
+ pt = f-> dog ;
+ for (s = s_min ; s <= s_max - 1 ; ++s) {
+ vl_sift_pix* src_a = vl_sift_get_octave (f, s ) ;
+ vl_sift_pix* src_b = vl_sift_get_octave (f, s + 1) ;
+ vl_sift_pix* end_a = src_a + w * h ;
+ while (src_a != end_a) {
+ *pt++ = *src_b++ - *src_a++ ;
+ }
+ }
+
+ /* -----------------------------------------------------------------
+ * Find local maxima of DoG
+ * -------------------------------------------------------------- */
+
+ /* start from dog [1,1,s_min+1] */
+ pt = dog + xo + yo + so ;
+
+ for(s = s_min + 1 ; s <= s_max - 2 ; ++s) {
+ for(y = 1 ; y < h - 1 ; ++y) {
+ for(x = 1 ; x < w - 1 ; ++x) {
+ v = *pt ;
+
+#define CHECK_NEIGHBORS(CMP,SGN) \
+ ( v CMP ## = SGN 0.8 * tp && \
+ v CMP *(pt + xo) && \
+ v CMP *(pt - xo) && \
+ v CMP *(pt + so) && \
+ v CMP *(pt - so) && \
+ v CMP *(pt + yo) && \
+ v CMP *(pt - yo) && \
+ \
+ v CMP *(pt + yo + xo) && \
+ v CMP *(pt + yo - xo) && \
+ v CMP *(pt - yo + xo) && \
+ v CMP *(pt - yo - xo) && \
+ \
+ v CMP *(pt + xo + so) && \
+ v CMP *(pt - xo + so) && \
+ v CMP *(pt + yo + so) && \
+ v CMP *(pt - yo + so) && \
+ v CMP *(pt + yo + xo + so) && \
+ v CMP *(pt + yo - xo + so) && \
+ v CMP *(pt - yo + xo + so) && \
+ v CMP *(pt - yo - xo + so) && \
+ \
+ v CMP *(pt + xo - so) && \
+ v CMP *(pt - xo - so) && \
+ v CMP *(pt + yo - so) && \
+ v CMP *(pt - yo - so) && \
+ v CMP *(pt + yo + xo - so) && \
+ v CMP *(pt + yo - xo - so) && \
+ v CMP *(pt - yo + xo - so) && \
+ v CMP *(pt - yo - xo - so) )
+
+ if (CHECK_NEIGHBORS(>,+) ||
+ CHECK_NEIGHBORS(<,-) ) {
+
+ /* make room for more keypoints */
+ if (f->nkeys >= f->keys_res) {
+ f->keys_res += 500 ;
+ if (f->keys) {
+ f->keys = vl_realloc (f->keys,
+ f->keys_res *
+ sizeof(VlSiftKeypoint)) ;
+ } else {
+ f->keys = vl_malloc (f->keys_res *
+ sizeof(VlSiftKeypoint)) ;
+ }
+ }
+
+ k = f->keys + (f->nkeys ++) ;
+
+ k-> ix = x ;
+ k-> iy = y ;
+ k-> is = s ;
+ }
+ pt += 1 ;
+ }
+ pt += 2 ;
+ }
+ pt += 2 * yo ;
+ }
+
+ /* -----------------------------------------------------------------
+ * Refine local maxima
+ * -------------------------------------------------------------- */
+
+ /* this pointer is used to write the keypoints back */
+ k = f->keys ;
+
+ for (i = 0 ; i < f->nkeys ; ++i) {
+
+ int x = f-> keys [i] .ix ;
+ int y = f-> keys [i] .iy ;
+ int s = f-> keys [i]. is ;
+
+ double Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
+ double A [3*3], b [3] ;
+
+ int dx = 0 ;
+ int dy = 0 ;
+
+ int iter, i, j ;
+
+ for (iter = 0 ; iter < 5 ; ++iter) {
+
+ x += dx ;
+ y += dy ;
+
+ pt = dog
+ + xo * x
+ + yo * y
+ + so * (s - s_min) ;
+
+ /** @brief Index GSS @internal */
+#define at(dx,dy,ds) (*( pt + (dx)*xo + (dy)*yo + (ds)*so))
+
+ /** @brief Index matrix A @internal */
+#define Aat(i,j) (A[(i)+(j)*3])
+
+ /* compute the gradient */
+ Dx = 0.5 * (at(+1,0,0) - at(-1,0,0)) ;
+ Dy = 0.5 * (at(0,+1,0) - at(0,-1,0));
+ Ds = 0.5 * (at(0,0,+1) - at(0,0,-1)) ;
+
+ /* compute the Hessian */
+ Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0 * at(0,0,0)) ;
+ Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0 * at(0,0,0)) ;
+ Dss = (at(0,0,+1) + at(0,0,-1) - 2.0 * at(0,0,0)) ;
+
+ Dxy = 0.25 * ( at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0) ) ;
+ Dxs = 0.25 * ( at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1) ) ;
+ Dys = 0.25 * ( at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1) ) ;
+
+ /* solve linear system ....................................... */
+ Aat(0,0) = Dxx ;
+ Aat(1,1) = Dyy ;
+ Aat(2,2) = Dss ;
+ Aat(0,1) = Aat(1,0) = Dxy ;
+ Aat(0,2) = Aat(2,0) = Dxs ;
+ Aat(1,2) = Aat(2,1) = Dys ;
+
+ b[0] = - Dx ;
+ b[1] = - Dy ;
+ b[2] = - Ds ;
+
+ /* Gauss elimination */
+ for(j = 0 ; j < 3 ; ++j) {
+ double maxa = 0 ;
+ double maxabsa = 0 ;
+ int maxi = -1 ;
+ double tmp ;
+
+ /* look for the maximally stable pivot */
+ for (i = j ; i < 3 ; ++i) {
+ double a = Aat (i,j) ;
+ double absa = vl_abs_d (a) ;
+ if (absa > maxabsa) {
+ maxa = a ;
+ maxabsa = absa ;
+ maxi = i ;
+ }
+ }
+
+ /* if singular give up */
+ if (maxabsa < 1e-10f) {
+ b[0] = 0 ;
+ b[1] = 0 ;
+ b[2] = 0 ;
+ break ;
+ }
+
+ i = maxi ;
+
+ /* swap j-th row with i-th row and normalize j-th row */
+ for(jj = j ; jj < 3 ; ++jj) {
+ tmp = Aat(i,jj) ; Aat(i,jj) = Aat(j,jj) ; Aat(j,jj) = tmp ;
+ Aat(j,jj) /= maxa ;
+ }
+ tmp = b[j] ; b[j] = b[i] ; b[i] = tmp ;
+ b[j] /= maxa ;
+
+ /* elimination */
+ for (ii = j+1 ; ii < 3 ; ++ii) {
+ double x = Aat(ii,j) ;
+ for (jj = j ; jj < 3 ; ++jj) {
+ Aat(ii,jj) -= x * Aat(j,jj) ;
+ }
+ b[ii] -= x * b[j] ;
+ }
+ }
+
+ /* backward substitution */
+ for (i = 2 ; i > 0 ; --i) {
+ double x = b[i] ;
+ for (ii = i-1 ; ii >= 0 ; --ii) {
+ b[ii] -= x * Aat(ii,i) ;
+ }
+ }
+
+ /* .......................................................... */
+ /* If the translation of the keypoint is big, move the keypoint
+ * and re-iterate the computation. Otherwise we are all set.
+ */
+
+ dx= ((b[0] > 0.6 && x < w - 2) ? 1 : 0)
+ + ((b[0] < -0.6 && x > 1 ) ? -1 : 0) ;
+
+ dy= ((b[1] > 0.6 && y < h - 2) ? 1 : 0)
+ + ((b[1] < -0.6 && y > 1 ) ? -1 : 0) ;
+
+ if (dx == 0 && dy == 0) break ;
+ }
+
+ /* check threshold and other conditions */
+ {
+ double val = at(0,0,0)
+ + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ;
+ double score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
+ double xn = x + b[0] ;
+ double yn = y + b[1] ;
+ double sn = s + b[2] ;
+
+ vl_bool good =
+ vl_abs_d (val) > tp &&
+ score < (te+1)*(te+1)/te &&
+ score >= 0 &&
+ vl_abs_d (b[0]) < 1.5 &&
+ vl_abs_d (b[1]) < 1.5 &&
+ vl_abs_d (b[2]) < 1.5 &&
+ xn >= 0 &&
+ xn <= w - 1 &&
+ yn >= 0 &&
+ yn <= h - 1 &&
+ sn >= s_min &&
+ sn <= s_max ;
+
+ if (good) {
+ k-> o = f->o_cur ;
+ k-> ix = x ;
+ k-> iy = y ;
+ k-> is = s ;
+ k-> s = sn ;
+ k-> x = xn * xper ;
+ k-> y = yn * xper ;
+ k-> sigma = f->sigma0 * pow (2.0, sn/f->S) * xper ;
+ ++ k ;
+ }
+
+ } /* done checking */
+ } /* next keypoint to refine */
+
+ /* update keypoint count */
+ f-> nkeys = k - f->keys ;
+}
+
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Update gradients to current GSS octave
+ **
+ ** @param f SIFT filter.
+ **
+ ** The function makes sure that the gradient buffer is up-to-date
+ ** with the current GSS data.
+ **
+ ** @remark The minimum octave size is 2x2xS.
+ **/
+
+static void
+update_gradient (VlSiftFilt *f)
+{
+ int s_min = f->s_min ;
+ int s_max = f->s_max ;
+ int w = vl_sift_get_octave_width (f) ;
+ int h = vl_sift_get_octave_height (f) ;
+ int const xo = 1 ;
+ int const yo = w ;
+ int const so = h * w ;
+ int y, s ;
+
+ if (f->grad_o == f->o_cur) return ;
+
+ for (s = s_min + 1 ;
+ s <= s_max - 2 ; ++ s) {
+
+ vl_sift_pix *src, *end, *grad, gx, gy ;
+
+#define SAVE_BACK \
+ *grad++ = vl_fast_sqrt_f (gx*gx + gy*gy) ; \
+ *grad++ = vl_mod_2pi_f (vl_fast_atan2_f (gy, gx) + 2*VL_PI) ; \
+ ++src ; \
+
+ src = vl_sift_get_octave (f,s) ;
+ grad = f->grad + 2 * so * (s - s_min -1) ;
+
+ /* first pixel of the first row */
+ gx = src[+xo] - src[0] ;
+ gy = src[+yo] - src[0] ;
+ SAVE_BACK ;
+
+ /* middle pixels of the first row */
+ end = (src - 1) + w - 1 ;
+ while (src < end) {
+ gx = 0.5 * (src[+xo] - src[-xo]) ;
+ gy = src[+yo] - src[0] ;
+ SAVE_BACK ;
+ }
+
+ /* last pixel of the first row */
+ gx = src[0] - src[-xo] ;
+ gy = src[+yo] - src[0] ;
+ SAVE_BACK ;
+
+ for (y = 1 ; y < h -1 ; ++y) {
+
+ /* first pixel of the middle rows */
+ gx = src[+xo] - src[0] ;
+ gy = 0.5 * (src[+yo] - src[-yo]) ;
+ SAVE_BACK ;
+
+ /* middle pixels of the middle rows */
+ end = (src - 1) + w - 1 ;
+ while (src < end) {
+ gx = 0.5 * (src[+xo] - src[-xo]) ;
+ gy = 0.5 * (src[+yo] - src[-yo]) ;
+ SAVE_BACK ;
+ }
+
+ /* last pixel of the middle row */
+ gx = src[0] - src[-xo] ;
+ gy = 0.5 * (src[+yo] - src[-yo]) ;
+ SAVE_BACK ;
+ }
+
+ /* first pixel of the last row */
+ gx = src[+xo] - src[0] ;
+ gy = src[ 0] - src[-yo] ;
+ SAVE_BACK ;
+
+ /* middle pixels of the last row */
+ end = (src - 1) + w - 1 ;
+ while (src < end) {
+ gx = 0.5 * (src[+xo] - src[-xo]) ;
+ gy = src[0] - src[-yo] ;
+ SAVE_BACK ;
+ }
+
+ /* last pixel of the last row */
+ gx = src[0] - src[-xo] ;
+ gy = src[0] - src[-yo] ;
+ SAVE_BACK ;
+ }
+ f->grad_o = f->o_cur ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Calculate the keypoint orientation(s)
+ **
+ ** @param f SIFT filter.
+ ** @param angles orientations (output).
+ ** @param k keypoint.
+ **
+ ** The function computes the orientation(s) of the keypoint @a k.
+ ** The function returns the number of orientations found (up to
+ ** four). The orientations themselves are written to the vector @a
+ ** angles.
+ **
+ ** @remark The function requires the keypoint octave @a k->o to be
+ ** equal to the filter current octave ::vl_sift_get_octave. If this
+ ** is not the case, the function returns zero orientations.
+ **
+ ** @remark The function requires the keypoint scale level @c k->s to
+ ** be in the range @c s_min+1 and @c s_max-2 (where usually @c
+ ** s_min=0 and @c s_max=S+2). If this is not the case, the function
+ ** returns zero orientations.
+ **
+ ** @return number of orientations found.
+ **/
+
+VL_EXPORT
+int
+vl_sift_calc_keypoint_orientations (VlSiftFilt *f,
+ double angles [4],
+ VlSiftKeypoint const *k)
+{
+ double const winf = 1.5 ;
+ double xper = pow (2.0, f->o_cur) ;
+
+ int w = f-> octave_width ;
+ int h = f-> octave_height ;
+ int const xo = 2 ; /* x-stride */
+ int const yo = 2 * w ; /* y-stride */
+ int const so = 2 * w * h ; /* s-stride */
+ double x = k-> x / xper ;
+ double y = k-> y / xper ;
+ double sigma = k-> sigma / xper ;
+
+ int xi = (int) (x + 0.5) ;
+ int yi = (int) (y + 0.5) ;
+ int si = k-> is ;
+
+ double const sigmaw = winf * sigma ;
+ int W = VL_MAX(floor (3.0 * sigmaw), 1) ;
+
+ int nangles= 0 ;
+
+ enum {nbins = 36} ;
+
+ double hist [nbins], maxh ;
+ vl_sift_pix const * pt ;
+ int xs, ys, iter, i ;
+
+ /* skip if the keypoint octave is not current */
+ if(k->o != f->o_cur)
+ return 0 ;
+
+ /* skip the keypoint if it is out of bounds */
+ if(xi < 0 ||
+ xi > w - 1 ||
+ yi < 0 ||
+ yi > h - 1 ||
+ si < f->s_min + 1 ||
+ si > f->s_max - 2 ) {
+ return 0 ;
+ }
+
+ /* make gradient up to date */
+ update_gradient (f) ;
+
+ /* clear histogram */
+ memset (hist, 0, sizeof(double) * nbins) ;
+
+ /* compute orientation histogram */
+ pt = f-> grad + xo*xi + yo*yi + so*(si - f->s_min - 1) ;
+
+#undef at
+#define at(dx,dy) (*(pt + xo * (dx) + yo * (dy)))
+
+ for(ys = VL_MAX (- W, - yi) ;
+ ys <= VL_MIN (+ W, h - 1 - yi) ; ++ys) {
+
+ for(xs = VL_MAX (- W, - xi) ;
+ xs <= VL_MIN (+ W, w - 1 - xi) ; ++xs) {
+
+
+ double dx = (double)(xi + xs) - x;
+ double dy = (double)(yi + ys) - y;
+ double r2 = dx*dx + dy*dy ;
+ double wgt, mod, ang, fbin ;
+
+ /* limit to a circular window */
+ if (r2 >= W*W + 0.6) continue ;
+
+ wgt = fast_expn (r2 / (2*sigmaw*sigmaw)) ;
+ mod = *(pt + xs*xo + ys*yo ) ;
+ ang = *(pt + xs*xo + ys*yo + 1) ;
+ fbin = nbins * ang / (2 * VL_PI) ;
+
+#if defined(VL_SIFT_BILINEAR_ORIENTATIONS)
+ {
+ int bin = vl_floor_d (fbin - 0.5) ;
+ double rbin = fbin - bin - 0.5 ;
+ hist [(bin + nbins) % nbins] += (1 - rbin) * mod * wgt ;
+ hist [(bin + 1 ) % nbins] += ( rbin) * mod * wgt ;
+ }
+#else
+ {
+ int bin = vl_floor_d (fbin) ;
+ bin = vl_floor_d (nbins * ang / (2*VL_PI)) ;
+ hist [(bin) % nbins] += mod * wgt ;
+ }
+#endif
+
+ } /* for xs */
+ } /* for ys */
+
+ /* smooth histogram */
+ for (iter = 0; iter < 6; iter ++) {
+ double prev = hist [nbins - 1] ;
+ double first = hist [0] ;
+ int i ;
+ for (i = 0; i < nbins - 1; i++) {
+ double newh = (prev + hist[i] + hist[(i+1) % nbins]) / 3.0;
+ prev = hist[i] ;
+ hist[i] = newh ;
+ }
+ hist[i] = (prev + hist[i] + first) / 3.0 ;
+ }
+
+ /* find the histogram maximum */
+ maxh = 0 ;
+ for (i = 0 ; i < nbins ; ++i)
+ maxh = VL_MAX (maxh, hist [i]) ;
+
+ /* find peaks within 80% from max */
+ nangles = 0 ;
+ for(i = 0 ; i < nbins ; ++i) {
+ double h0 = hist [i] ;
+ double hm = hist [(i - 1 + nbins) % nbins] ;
+ double hp = hist [(i + 1 + nbins) % nbins] ;
+
+ /* is this a peak? */
+ if (h0 > 0.8*maxh && h0 > hm && h0 > hp) {
+
+ /* quadratic interpolation */
+ double di = - 0.5 * (hp - hm) / (hp + hm - 2 * h0) ;
+ double th = 2 * VL_PI * (i + di + 0.5) / nbins ;
+ angles [ nangles++ ] = th ;
+ if( nangles == 4 )
+ goto enough_angles ;
+ }
+ }
+ enough_angles:
+ return nangles ;
+}
+
+
+/** ------------------------------------------------------------------
+ ** @internal
+ ** @brief Normalizes in norm L_2 a descriptor
+ ** @param begin begin of histogram.
+ ** @param end end of histogram.
+ **/
+
+VL_INLINE vl_sift_pix
+normalize_histogram
+(vl_sift_pix *begin, vl_sift_pix *end)
+{
+ vl_sift_pix* iter ;
+ vl_sift_pix norm = 0.0 ;
+
+ for (iter = begin ; iter != end ; ++ iter)
+ norm += (*iter) * (*iter) ;
+
+ norm = vl_fast_sqrt_f (norm) + VL_EPSILON_F ;
+
+ for (iter = begin; iter != end ; ++ iter)
+ *iter /= norm ;
+
+ return norm;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Run the SIFT descriptor on raw data
+ **
+ ** @param f SIFT filter.
+ ** @param grad image gradients.
+ ** @param descr SIFT descriptor (output).
+ ** @param width image width.
+ ** @param height image height.
+ ** @param x keypoint x coordinate.
+ ** @param y keypoint y coordinate.
+ ** @param sigma keypoint scale.
+ ** @param angle0 keypoint orientation.
+ **
+ ** The function runs the SIFT descriptor on raw data. Here @a image
+ ** is a 2 x @a width x @a height array (by convention, the memory
+ ** layout is a s such the first index is the fastest varying
+ ** one). The first @a width x @a height layer of the array contains
+ ** the gradient magnitude and the second the gradient angle (in
+ ** radians, between 0 and @f$ 2\pi @f$). @a x, @a y and @a sigma give
+ ** the keypoint center and scale respectively.
+ **
+ ** In order to be equivalent to a standard SIFT descriptor the image
+ ** gradient must be computed at a smoothing level equal to the scale
+ ** of the keypoint. In practice, the actual SIFT algorithm makes the
+ ** following additional approximation, which influence the result:
+ **
+ ** - Scale is discretized in @c S levels.
+ ** - The image is downsampled once for each octave (if you do this,
+ ** the parameters @a x, @a y and @a sigma must be
+ ** scaled too).
+ **/
+
+VL_EXPORT
+void
+vl_sift_calc_raw_descriptor (VlSiftFilt const *f,
+ vl_sift_pix const* grad,
+ vl_sift_pix *descr,
+ int width, int height,
+ double x, double y,
+ double sigma,
+ double angle0)
+{
+ double const magnif = f-> magnif ;
+
+ int w = width ;
+ int h = height ;
+ int const xo = 2 ; /* x-stride */
+ int const yo = 2 * w ; /* y-stride */
+
+ int xi = (int) (x + 0.5) ;
+ int yi = (int) (y + 0.5) ;
+
+ double const st0 = sin (angle0) ;
+ double const ct0 = cos (angle0) ;
+ double const SBP = magnif * sigma + VL_EPSILON_D ;
+ int const W = floor
+ (sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5) ;
+
+ int const binto = 1 ; /* bin theta-stride */
+ int const binyo = NBO * NBP ; /* bin y-stride */
+ int const binxo = NBO ; /* bin x-stride */
+
+ int bin, dxi, dyi ;
+ vl_sift_pix const *pt ;
+ vl_sift_pix *dpt ;
+
+ /* check bounds */
+ if(xi < 0 ||
+ xi >= w ||
+ yi < 0 ||
+ yi >= h - 1 )
+ return ;
+
+ /* clear descriptor */
+ memset (descr, 0, sizeof(vl_sift_pix) * NBO*NBP*NBP) ;
+
+ /* Center the scale space and the descriptor on the current keypoint.
+ * Note that dpt is pointing to the bin of center (SBP/2,SBP/2,0).
+ */
+ pt = grad + xi*xo + yi*yo ;
+ dpt = descr + (NBP/2) * binyo + (NBP/2) * binxo ;
+
+#undef atd
+#define atd(dbinx,dbiny,dbint) *(dpt + (dbint)*binto + (dbiny)*binyo + (dbinx)*binxo)
+
+ /*
+ * Process pixels in the intersection of the image rectangle
+ * (1,1)-(M-1,N-1) and the keypoint bounding box.
+ */
+ for(dyi = VL_MAX(- W, - yi ) ;
+ dyi <= VL_MIN(+ W, h - yi -1) ; ++ dyi) {
+
+ for(dxi = VL_MAX(- W, - xi ) ;
+ dxi <= VL_MIN(+ W, w - xi -1) ; ++ dxi) {
+
+ /* retrieve */
+ vl_sift_pix mod = *( pt + dxi*xo + dyi*yo + 0 ) ;
+ vl_sift_pix angle = *( pt + dxi*xo + dyi*yo + 1 ) ;
+ vl_sift_pix theta = vl_mod_2pi_f (angle - angle0) ;
+
+ /* fractional displacement */
+ vl_sift_pix dx = xi + dxi - x;
+ vl_sift_pix dy = yi + dyi - y;
+
+ /* get the displacement normalized w.r.t. the keypoint
+ orientation and extension */
+ vl_sift_pix nx = ( ct0 * dx + st0 * dy) / SBP ;
+ vl_sift_pix ny = (-st0 * dx + ct0 * dy) / SBP ;
+ vl_sift_pix nt = NBO * theta / (2 * VL_PI) ;
+
+ /* Get the Gaussian weight of the sample. The Gaussian window
+ * has a standard deviation equal to NBP/2. Note that dx and dy
+ * are in the normalized frame, so that -NBP/2 <= dx <=
+ * NBP/2. */
+ vl_sift_pix const wsigma = f->windowSize ;
+ vl_sift_pix win = fast_expn
+ ((nx*nx + ny*ny)/(2.0 * wsigma * wsigma)) ;
+
+ /* The sample will be distributed in 8 adjacent bins.
+ We start from the ``lower-left'' bin. */
+ int binx = vl_floor_f (nx - 0.5) ;
+ int biny = vl_floor_f (ny - 0.5) ;
+ int bint = vl_floor_f (nt) ;
+ vl_sift_pix rbinx = nx - (binx + 0.5) ;
+ vl_sift_pix rbiny = ny - (biny + 0.5) ;
+ vl_sift_pix rbint = nt - bint ;
+ int dbinx ;
+ int dbiny ;
+ int dbint ;
+
+ /* Distribute the current sample into the 8 adjacent bins*/
+ for(dbinx = 0 ; dbinx < 2 ; ++dbinx) {
+ for(dbiny = 0 ; dbiny < 2 ; ++dbiny) {
+ for(dbint = 0 ; dbint < 2 ; ++dbint) {
+
+ if (binx + dbinx >= - (NBP/2) &&
+ binx + dbinx < (NBP/2) &&
+ biny + dbiny >= - (NBP/2) &&
+ biny + dbiny < (NBP/2) ) {
+ vl_sift_pix weight = win
+ * mod
+ * vl_abs_f (1 - dbinx - rbinx)
+ * vl_abs_f (1 - dbiny - rbiny)
+ * vl_abs_f (1 - dbint - rbint) ;
+
+ atd(binx+dbinx, biny+dbiny, (bint+dbint) % NBO) += weight ;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* Standard SIFT descriptors are normalized, truncated and normalized again */
+ if(1) {
+
+ /* normalize L2 norm */
+ vl_sift_pix norm = normalize_histogram (descr, descr + NBO*NBP*NBP) ;
+
+ /*
+ Set the descriptor to zero if it is lower than our
+ norm_threshold. We divide by the number of samples in the
+ descriptor region becasue the Gaussian window used in the
+ calculation fo the descriptor is not normalized.
+ */
+ int numSamples =
+ (VL_MIN(W, w - xi -1) - VL_MAX(-W, - xi) + 1) *
+ (VL_MIN(W, h - yi -1) - VL_MAX(-W, - yi) + 1) ;
+
+ if(f-> norm_thresh && norm < f-> norm_thresh * numSamples) {
+ for(bin = 0; bin < NBO*NBP*NBP ; ++ bin)
+ descr [bin] = 0;
+ }
+ else {
+ /* truncate at 0.2. */
+ for(bin = 0; bin < NBO*NBP*NBP ; ++ bin) {
+ if (descr [bin] > 0.2) descr [bin] = 0.2;
+ }
+
+ /* normalize again. */
+ normalize_histogram (descr, descr + NBO*NBP*NBP) ;
+ }
+ }
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Compute the descriptor of a keypoint
+ **
+ ** @param f SIFT filter.
+ ** @param descr SIFT descriptor (output)
+ ** @param k keypoint.
+ ** @param angle0 keypoint direction.
+ **
+ ** The function computes the SIFT descriptor of the keypoint @a k of
+ ** orientation @a angle0. The function fills the buffer @a descr
+ ** which must be large enough to hold the descriptor.
+ **
+ ** The function assumes that the keypoint is on the current octave.
+ ** If not, it does not do anything.
+ **/
+
+VL_EXPORT
+void
+vl_sift_calc_keypoint_descriptor (VlSiftFilt *f,
+ vl_sift_pix *descr,
+ VlSiftKeypoint const* k,
+ double angle0)
+{
+ /*
+ The SIFT descriptor is a three dimensional histogram of the
+ position and orientation of the gradient. There are NBP bins for
+ each spatial dimension and NBO bins for the orientation dimension,
+ for a total of NBP x NBP x NBO bins.
+
+ The support of each spatial bin has an extension of SBP = 3sigma
+ pixels, where sigma is the scale of the keypoint. Thus all the
+ bins together have a support SBP x NBP pixels wide. Since
+ weighting and interpolation of pixel is used, the support extends
+ by another half bin. Therefore, the support is a square window of
+ SBP x (NBP + 1) pixels. Finally, since the patch can be
+ arbitrarily rotated, we need to consider a window 2W += sqrt(2) x
+ SBP x (NBP + 1) pixels wide.
+ */
+
+ double const magnif = f-> magnif ;
+
+ double xper = pow (2.0, f->o_cur) ;
+
+ int w = f-> octave_width ;
+ int h = f-> octave_height ;
+ int const xo = 2 ; /* x-stride */
+ int const yo = 2 * w ; /* y-stride */
+ int const so = 2 * w * h ; /* s-stride */
+ double x = k-> x / xper ;
+ double y = k-> y / xper ;
+ double sigma = k-> sigma / xper ;
+
+ int xi = (int) (x + 0.5) ;
+ int yi = (int) (y + 0.5) ;
+ int si = k-> is ;
+
+ double const st0 = sin (angle0) ;
+ double const ct0 = cos (angle0) ;
+ double const SBP = magnif * sigma ;
+ int const W = floor
+ (sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5) ;
+
+ int const binto = 1 ; /* bin theta-stride */
+ int const binyo = NBO * NBP ; /* bin y-stride */
+ int const binxo = NBO ; /* bin x-stride */
+
+ int bin, dxi, dyi ;
+ vl_sift_pix const *pt ;
+ vl_sift_pix *dpt ;
+
+ /* check bounds */
+ if(k->o != f->o_cur ||
+ xi < 0 ||
+ xi >= w ||
+ yi < 0 ||
+ yi >= h - 1 ||
+ si < f->s_min + 1 ||
+ si > f->s_max - 2 )
+ return ;
+
+ /* synchronize gradient buffer */
+ update_gradient (f) ;
+
+ /* VL_PRINTF("W = %d ; magnif = %g ; SBP = %g\n", W,magnif,SBP) ; */
+
+ /* clear descriptor */
+ memset (descr, 0, sizeof(vl_sift_pix) * NBO*NBP*NBP) ;
+
+ /* Center the scale space and the descriptor on the current keypoint.
+ * Note that dpt is pointing to the bin of center (SBP/2,SBP/2,0).
+ */
+ pt = f->grad + xi*xo + yi*yo + (si - f->s_min - 1)*so ;
+ dpt = descr + (NBP/2) * binyo + (NBP/2) * binxo ;
+
+#undef atd
+#define atd(dbinx,dbiny,dbint) *(dpt + (dbint)*binto + (dbiny)*binyo + (dbinx)*binxo)
+
+ /*
+ * Process pixels in the intersection of the image rectangle
+ * (1,1)-(M-1,N-1) and the keypoint bounding box.
+ */
+ for(dyi = VL_MAX (- W, 1 - yi ) ;
+ dyi <= VL_MIN (+ W, h - yi - 2) ; ++ dyi) {
+
+ for(dxi = VL_MAX (- W, 1 - xi ) ;
+ dxi <= VL_MIN (+ W, w - xi - 2) ; ++ dxi) {
+
+ /* retrieve */
+ vl_sift_pix mod = *( pt + dxi*xo + dyi*yo + 0 ) ;
+ vl_sift_pix angle = *( pt + dxi*xo + dyi*yo + 1 ) ;
+ vl_sift_pix theta = vl_mod_2pi_f (angle - angle0) ;
+
+ /* fractional displacement */
+ vl_sift_pix dx = xi + dxi - x;
+ vl_sift_pix dy = yi + dyi - y;
+
+ /* get the displacement normalized w.r.t. the keypoint
+ orientation and extension */
+ vl_sift_pix nx = ( ct0 * dx + st0 * dy) / SBP ;
+ vl_sift_pix ny = (-st0 * dx + ct0 * dy) / SBP ;
+ vl_sift_pix nt = NBO * theta / (2 * VL_PI) ;
+
+ /* Get the Gaussian weight of the sample. The Gaussian window
+ * has a standard deviation equal to NBP/2. Note that dx and dy
+ * are in the normalized frame, so that -NBP/2 <= dx <=
+ * NBP/2. */
+ vl_sift_pix const wsigma = f->windowSize ;
+ vl_sift_pix win = fast_expn
+ ((nx*nx + ny*ny)/(2.0 * wsigma * wsigma)) ;
+
+ /* The sample will be distributed in 8 adjacent bins.
+ We start from the ``lower-left'' bin. */
+ int binx = vl_floor_f (nx - 0.5) ;
+ int biny = vl_floor_f (ny - 0.5) ;
+ int bint = vl_floor_f (nt) ;
+ vl_sift_pix rbinx = nx - (binx + 0.5) ;
+ vl_sift_pix rbiny = ny - (biny + 0.5) ;
+ vl_sift_pix rbint = nt - bint ;
+ int dbinx ;
+ int dbiny ;
+ int dbint ;
+
+ /* Distribute the current sample into the 8 adjacent bins*/
+ for(dbinx = 0 ; dbinx < 2 ; ++dbinx) {
+ for(dbiny = 0 ; dbiny < 2 ; ++dbiny) {
+ for(dbint = 0 ; dbint < 2 ; ++dbint) {
+
+ if (binx + dbinx >= - (NBP/2) &&
+ binx + dbinx < (NBP/2) &&
+ biny + dbiny >= - (NBP/2) &&
+ biny + dbiny < (NBP/2) ) {
+ vl_sift_pix weight = win
+ * mod
+ * vl_abs_f (1 - dbinx - rbinx)
+ * vl_abs_f (1 - dbiny - rbiny)
+ * vl_abs_f (1 - dbint - rbint) ;
+
+ atd(binx+dbinx, biny+dbiny, (bint+dbint) % NBO) += weight ;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* Standard SIFT descriptors are normalized, truncated and normalized again */
+ if(1) {
+
+ /* Normalize the histogram to L2 unit length. */
+ vl_sift_pix norm = normalize_histogram (descr, descr + NBO*NBP*NBP) ;
+
+ /* Set the descriptor to zero if it is lower than our norm_threshold */
+ if(f-> norm_thresh && norm < f-> norm_thresh) {
+ for(bin = 0; bin < NBO*NBP*NBP ; ++ bin)
+ descr [bin] = 0;
+ }
+ else {
+
+ /* Truncate at 0.2. */
+ for(bin = 0; bin < NBO*NBP*NBP ; ++ bin) {
+ if (descr [bin] > 0.2) descr [bin] = 0.2;
+ }
+
+ /* Normalize again. */
+ normalize_histogram (descr, descr + NBO*NBP*NBP) ;
+ }
+ }
+
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Initialize a keypoint from its position and scale
+ **
+ ** @param f SIFT filter.
+ ** @param k SIFT keypoint (output).
+ ** @param x x coordinate of the center.
+ ** @param y y coordinate of the center.
+ ** @param sigma scale.
+ **
+ ** The function initializes the structure @a k from the location @a x
+ ** and @a y and scale @a sigma of the keypoint.
+ **/
+
+VL_EXPORT
+void
+vl_sift_keypoint_init (VlSiftFilt const *f,
+ VlSiftKeypoint *k,
+ double x,
+ double y,
+ double sigma)
+{
+
+ /*
+ The formula linking the keypoint scale sigma to the octave and
+ scale index is
+
+ (1) sigma(o,s) = sigma0 2^(o+s/S)
+
+ for which
+
+ (2) o + s/S = log2 sigma/sigma0 == phi.
+
+ In addition to the scale index s (which can be fractional due to
+ scale interpolation) a keypoint has an integer scale index is too
+ (which is the index of the scale level where it was detected in
+ the DoG scale space). We have the constraints:
+
+ - o and is are integer
+
+ - is is in the range [smin+1, smax-2 ]
+
+ - o is in the range [omin, omin+O-1]
+
+ - is = rand(s) most of the times (but not always, due to the way s
+ is obtained by quadratic interpolation of the DoG scale space).
+
+ Depending on the values of smin and smax, often (2) has multiple
+ solutions is,o that satisfy all constraints. In this case we
+ choose the one with biggest index o (this saves a bit of
+ computation).
+
+ DETERMINING THE OCTAVE INDEX O
+
+ From (2) we have o = phi - s/S and we want to pick the biggest
+ possible index o in the feasible range. This corresponds to
+ selecting the smallest possible index s. We write s = is + ds
+ where in most cases |ds|<.5 (but in general |ds|<1). So we have
+
+ o = phi - s/S, s = is + ds , |ds| < .5 (or |ds| < 1).
+
+ Since is is in the range [smin+1,smax-2], s is in the range
+ [smin+.5,smax-1.5] (or [smin,smax-1]), the number o is an integer
+ in the range phi+[-smax+1.5,-smin-.5] (or
+ phi+[-smax+1,-smin]). Thus the maximum value of o is obtained for
+ o = floor(phi-smin-.5) (or o = floor(phi-smin)).
+
+ Finally o is clamped to make sure it is contained in the feasible
+ range.
+
+ DETERMINING THE SCALE INDEXES S AND IS
+
+ Given o we can derive is by writing (2) as
+
+ s = is + ds = S(phi - o).
+
+ We then take is = round(s) and clamp its value to be in the
+ feasible range.
+ */
+
+ int o, ix, iy, is ;
+ double s, phi, xper ;
+
+ phi = log2 (sigma / f->sigma0) ;
+ o = vl_floor_d (phi - ((double) f->s_min + 0.5) / f->S) ;
+ o = VL_MIN (o, f->o_min + f->O - 1) ;
+ o = VL_MAX (o, f->o_min ) ;
+ s = f->S * (phi - o) ;
+
+ is = (int)(s + 0.5) ;
+ is = VL_MIN(is, f->s_max - 2) ;
+ is = VL_MAX(is, f->s_min + 1) ;
+
+ xper = pow (2.0, o) ;
+ ix = (int)(x / xper + 0.5) ;
+ iy = (int)(y / xper + 0.5) ;
+
+ k -> o = o ;
+
+ k -> ix = ix ;
+ k -> iy = iy ;
+ k -> is = is ;
+
+ k -> x = x ;
+ k -> y = y ;
+ k -> s = s ;
+
+ k->sigma = sigma ;
+
+ /*
+ VL_PRINTF ("k.ix = %d\n", ix) ;
+ VL_PRINTF ("k.iy = %d\n", iy) ;
+ VL_PRINTF ("k.is = %d\n", is) ;
+ VL_PRINTF ("k.o = %d\n", o ) ;
+ VL_PRINTF ("k.s = %g\n", s ) ;
+ VL_PRINTF ("k.x = %g\n", x ) ;
+ VL_PRINTF ("k.y = %g\n", y ) ;
+ VL_PRINTF ("k.sigma = %g\n", sigma) ;
+ */
+}
diff --git a/vl/sift.h b/vl/sift.h
new file mode 100644
index 0000000..1db4668
--- /dev/null
+++ b/vl/sift.h
@@ -0,0 +1,410 @@
+/** @file sift.h
+ ** @brief Scale Invariant Feature Transform (SIFT)
+ ** @author Andrea Vedaldi
+ **/
+
+/* AUTORIGHTS
+Copyright 2007 (c) Andrea Vedaldi and Brian Fulkerson
+
+This file is part of VLFeat, available in the terms of the GNU
+General Public License version 2.
+*/
+
+#ifndef VL_SIFT_H
+#define VL_SIFT_H
+
+#include <stdio.h>
+#include "generic.h"
+
+/** @brief SIFT filter pixel type */
+typedef float vl_sift_pix ;
+
+/** ------------------------------------------------------------------
+ ** @brief SIFT filter keypoint
+ **
+ ** This structure represent a keypoint as extracted by the SIFT
+ ** filter ::VLSiftFilt.
+ **/
+
+typedef struct _VlSiftKeypoint
+{
+ int o ; /**< o coordinate (octave). */
+
+ int ix ; /**< Integer unnormalized x coordinate. */
+ int iy ; /**< Integer unnormalized y coordinate. */
+ int is ; /**< Integer s coordinate. */
+
+ float x ; /**< x coordinate. */
+ float y ; /**< u coordinate. */
+ float s ; /**< x coordinate. */
+ float sigma ; /**< scale. */
+} VlSiftKeypoint ;
+
+/** ------------------------------------------------------------------
+ ** @brief SIFT filter
+ **
+ ** This filter implements the SIFT detector and descriptor.
+ **/
+
+typedef struct _VlSiftFilt
+{
+ double sigman ; /**< nominal image smoothing. */
+ double sigma0 ; /**< smoothing of pyramid base. */
+ double sigmak ; /**< k-smoothing */
+ double dsigma0 ; /**< delta-smoothing. */
+
+ int width ; /**< image width. */
+ int height ; /**< image height. */
+ int O ; /**< number of octaves. */
+ int S ; /**< number of levels per octave. */
+ int o_min ; /**< minimum octave index. */
+ int s_min ; /**< minimum level index. */
+ int s_max ; /**< maximum level index. */
+ int o_cur ; /**< current octave. */
+
+ vl_sift_pix *temp ; /**< temporary pixel buffer. */
+ vl_sift_pix *octave ; /**< current GSS data. */
+ vl_sift_pix *dog ; /**< current DoG data. */
+ int octave_width ; /**< current octave width. */
+ int octave_height ; /**< current octave height. */
+
+ VlSiftKeypoint* keys ;/**< detected keypoints. */
+ int nkeys ; /**< number of detected keypoints. */
+ int keys_res ; /**< size of the keys buffer. */
+
+ double peak_thresh ; /**< peak threshold. */
+ double edge_thresh ; /**< edge threshold. */
+ double norm_thresh ; /**< norm threshold. */
+ double magnif ; /**< magnification factor. */
+ double windowSize ; /**< size of Gaussian window (in spatial bins) */
+
+ vl_sift_pix *grad ; /**< GSS gradient data. */
+ int grad_o ; /**< GSS gradient data octave. */
+
+} VlSiftFilt ;
+
+/** @name Create and destroy
+ ** @{
+ **/
+VL_EXPORT
+VlSiftFilt* vl_sift_new (int width, int height,
+ int noctaves, int nlevels,
+ int o_min) ;
+VL_EXPORT
+void vl_sift_delete (VlSiftFilt *f) ;
+/** @} */
+
+/** @name Process data
+ ** @{
+ **/
+
+VL_EXPORT
+int vl_sift_process_first_octave (VlSiftFilt *f,
+ vl_sift_pix const *im) ;
+
+VL_EXPORT
+int vl_sift_process_next_octave (VlSiftFilt *f) ;
+
+VL_EXPORT
+void vl_sift_detect (VlSiftFilt *f) ;
+
+VL_EXPORT
+int vl_sift_calc_keypoint_orientations (VlSiftFilt *f,
+ double angles [4],
+ VlSiftKeypoint const*k);
+VL_EXPORT
+void vl_sift_calc_keypoint_descriptor (VlSiftFilt *f,
+ vl_sift_pix *descr,
+ VlSiftKeypoint const* k,
+ double angle) ;
+
+VL_EXPORT
+void vl_sift_calc_raw_descriptor (VlSiftFilt const *f,
+ vl_sift_pix const* image,
+ vl_sift_pix *descr,
+ int widht, int height,
+ double x, double y,
+ double s, double angle0) ;
+
+VL_EXPORT
+void vl_sift_keypoint_init (VlSiftFilt const *f,
+ VlSiftKeypoint *k,
+ double x,
+ double y,
+ double sigma) ;
+/** @} */
+
+/** @name Retrieve data and parameters
+ ** @{
+ **/
+VL_INLINE int vl_sift_get_octave_index (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_noctaves (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_octave_first (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_octave_width (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_octave_height (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_nlevels (VlSiftFilt const *f) ;
+VL_INLINE int vl_sift_get_nkeypoints (VlSiftFilt const *f) ;
+VL_INLINE double vl_sift_get_peak_thresh (VlSiftFilt const *f) ;
+VL_INLINE double vl_sift_get_edge_thresh (VlSiftFilt const *f) ;
+VL_INLINE double vl_sift_get_norm_thresh (VlSiftFilt const *f) ;
+VL_INLINE double vl_sift_get_magnif (VlSiftFilt const *f) ;
+VL_INLINE double vl_sift_get_window_size (VlSiftFilt const *f) ;
+
+VL_INLINE vl_sift_pix *vl_sift_get_octave (VlSiftFilt const *f, int s) ;
+VL_INLINE VlSiftKeypoint const *vl_sift_get_keypoints (VlSiftFilt const *f) ;
+/** @} */
+
+/** @name Set parameters
+ ** @{
+ **/
+VL_INLINE void vl_sift_set_peak_thresh (VlSiftFilt *f, double t) ;
+VL_INLINE void vl_sift_set_edge_thresh (VlSiftFilt *f, double t) ;
+VL_INLINE void vl_sift_set_norm_thresh (VlSiftFilt *f, double t) ;
+VL_INLINE void vl_sift_set_magnif (VlSiftFilt *f, double m) ;
+VL_INLINE void vl_sift_set_window_size (VlSiftFilt *f, double m) ;
+/** @} */
+
+/* -------------------------------------------------------------------
+ * Inline functions implementation
+ * ---------------------------------------------------------------- */
+
+/** ------------------------------------------------------------------
+ ** @brief Get current octave index.
+ ** @param f SIFT filter.
+ ** @return index of the current octave.
+ **/
+
+VL_INLINE int
+vl_sift_get_octave_index (VlSiftFilt const *f)
+{
+ return f-> o_cur ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get number of octaves.
+ ** @param f SIFT filter.
+ ** @return number of octaves.
+ **/
+
+VL_INLINE int
+vl_sift_get_noctaves (VlSiftFilt const *f)
+{
+ return f-> O ;
+}
+
+/**-------------------------------------------------------------------
+ ** @brief Get first octave.
+ ** @param f SIFT filter.
+ ** @return index of the first octave.
+ **/
+
+VL_INLINE int
+vl_sift_get_octave_first (VlSiftFilt const *f)
+{
+ return f-> o_min ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get current octave width
+ ** @param f SIFT filter.
+ ** @return current octave width.
+ **/
+
+VL_INLINE int
+vl_sift_get_octave_width (VlSiftFilt const *f)
+{
+ return f-> octave_width ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get current octave height
+ ** @param f SIFT filter.
+ ** @return current octave height.
+ **/
+
+VL_INLINE int
+vl_sift_get_octave_height (VlSiftFilt const *f)
+{
+ return f-> octave_height ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get current octave data
+ ** @param f SIFT filter.
+ ** @param s level index.
+ **
+ ** The level index @a s ranges in the interval <tt>s_min = -1</tt>
+ ** and <tt> s_max = S + 2</tt>, where @c S is the number of levels
+ ** per octave.
+ **
+ ** @return pointer to the octave data for level @a s.
+ **/
+
+VL_INLINE vl_sift_pix *
+vl_sift_get_octave (VlSiftFilt const *f, int s)
+{
+ int w = vl_sift_get_octave_width (f) ;
+ int h = vl_sift_get_octave_height (f) ;
+ return f->octave + w * h * (s - f->s_min) ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get number of levels per octave
+ ** @param f SIFT filter.
+ ** @return number of leves per octave.
+ **/
+
+VL_INLINE int
+vl_sift_get_nlevels (VlSiftFilt const *f)
+{
+ return f-> S ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get number of keypoints.
+ ** @param f SIFT filter.
+ ** @return number of keypoints.
+ **/
+
+VL_INLINE int
+vl_sift_get_nkeypoints (VlSiftFilt const *f)
+{
+ return f-> nkeys ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get keypoints.
+ ** @param f SIFT filter.
+ ** @return pointer to the keypoints list.
+ **/
+
+VL_INLINE VlSiftKeypoint const *
+vl_sift_get_keypoints (VlSiftFilt const *f)
+{
+ return f-> keys ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get peaks treashold
+ ** @param f SIFT filter.
+ ** @return threshold ;
+ **/
+
+VL_INLINE double
+vl_sift_get_peak_thresh (VlSiftFilt const *f)
+{
+ return f -> peak_thresh ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get edges threshold
+ ** @param f SIFT filter.
+ ** @return threshold.
+ **/
+
+VL_INLINE double
+vl_sift_get_edge_thresh (VlSiftFilt const *f)
+{
+ return f -> edge_thresh ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get norm threshold
+ ** @param f SIFT filter.
+ ** @return threshold.
+ **/
+
+VL_INLINE double
+vl_sift_get_norm_thresh (VlSiftFilt const *f)
+{
+ return f -> norm_thresh ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get the magnification factor
+ ** @param f SIFT filter.
+ ** @return magnification factor.
+ **/
+
+VL_INLINE double
+vl_sift_get_magnif (VlSiftFilt const *f)
+{
+ return f -> magnif ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Get the Gaussian window size.
+ ** @param f SIFT filter.
+ ** @return standard deviation of the Gaussian window (in spatial bin units).
+ **/
+
+VL_INLINE double
+vl_sift_get_window_size (VlSiftFilt const *f)
+{
+ return f -> windowSize ;
+}
+
+
+
+/** ------------------------------------------------------------------
+ ** @brief Set peaks threshold
+ ** @param f SIFT filter.
+ ** @param t threshold.
+ **/
+
+VL_INLINE void
+vl_sift_set_peak_thresh (VlSiftFilt *f, double t)
+{
+ f -> peak_thresh = t ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Set edges threshold
+ ** @param f SIFT filter.
+ ** @param t threshold.
+ **/
+
+VL_INLINE void
+vl_sift_set_edge_thresh (VlSiftFilt *f, double t)
+{
+ f -> edge_thresh = t ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Set norm threshold
+ ** @param f SIFT filter.
+ ** @param t threshold.
+ **/
+
+VL_INLINE void
+vl_sift_set_norm_thresh (VlSiftFilt *f, double t)
+{
+ f -> norm_thresh = t ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Set the magnification factor
+ ** @param f SIFT filter.
+ ** @param m magnification factor.
+ **/
+
+VL_INLINE void
+vl_sift_set_magnif (VlSiftFilt *f, double m)
+{
+ f -> magnif = m ;
+}
+
+/** ------------------------------------------------------------------
+ ** @brief Set the Gaussian window size
+ ** @param f SIFT filter.
+ ** @param x Gaussian window size (in spatial bin units).
+ **/
+
+VL_INLINE void
+vl_sift_set_window_size (VlSiftFilt *f, double x)
+{
+ f -> windowSize = x ;
+}
+
+/* VL_SIFT_H */
+#endif