summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMichael M <myrhev@Myrhev-laptop.(none)>2009-11-19 13:45:27 +0100
committerMichael M <myrhev@Myrhev-laptop.(none)>2009-11-19 13:45:27 +0100
commita74794f28106dc16c2158c3d20ad3e6aa51bd4ac (patch)
treeef75992205ebf6448658917920fe438f0810d7c9
parentc503c0dd53d74c8807be0a75e0c0e4ccec77a935 (diff)
downloadpanor-a74794f28106dc16c2158c3d20ad3e6aa51bd4ac.tar.gz
panor-a74794f28106dc16c2158c3d20ad3e6aa51bd4ac.tar.bz2
Missing files
-rw-r--r--SIFT.cpp53
-rw-r--r--framework.cpp127
-rw-r--r--framework.h32
3 files changed, 212 insertions, 0 deletions
diff --git a/SIFT.cpp b/SIFT.cpp
new file mode 100644
index 0000000..7d0488a
--- /dev/null
+++ b/SIFT.cpp
@@ -0,0 +1,53 @@
+#include <iostream>
+#include <list>
+
+#include "Imagine/Features.h"
+
+using namespace std;
+
+extern "C" {
+#include "vl/sift.h"
+}
+
+namespace Imagine {
+
+ Array<SIFT> SIFTDetector::run(const Image<byte>& I) const {
+ int w = I.width(), h = I.height();
+ Image<float,2> If(I);
+
+ VlSiftFilt* filt = vl_sift_new(w, h, numOctaves, numScales, firstOctave);
+ if (edgeThresh >= 0)
+ vl_sift_set_edge_thresh(filt, edgeThresh);
+ if (peakThresh >= 0)
+ vl_sift_set_peak_thresh(filt, 255 * peakThresh / numScales);
+
+ list<SIFT> L;
+ vl_sift_process_first_octave(filt, If.data());
+ while (true) {
+ vl_sift_detect(filt);
+
+ VlSiftKeypoint const *keys = vl_sift_get_keypoints(filt) ;
+ int nkeys = vl_sift_get_nkeypoints(filt) ;
+ for(int i = 0; i < nkeys; ++i) {
+ double angles[4];
+ int nangles = vl_sift_calc_keypoint_orientations(filt, angles, keys+i) ;
+ for (int q = 0; q < nangles; ++q) {
+ vl_sift_pix descr[128] ;
+ vl_sift_calc_keypoint_descriptor(filt, descr, keys+i, angles[q]) ;
+ SIFT fp;
+ fp.pos = fPoint2(keys[i].x, keys[i].y);
+ fp.scale = keys[i].sigma;
+ fp.angle = float(angles[q]);
+ for (int k = 0; k < 128; k++)
+ fp.desc[k] = byte(512*descr[k]);
+ L.push_back(fp);
+ }
+ }
+ if (vl_sift_process_next_octave(filt))
+ break; // Last octave
+ }
+ vl_sift_delete(filt);
+ return(Array<SIFT>(L));
+ }
+
+}
diff --git a/framework.cpp b/framework.cpp
new file mode 100644
index 0000000..9c28afc
--- /dev/null
+++ b/framework.cpp
@@ -0,0 +1,127 @@
+#include "framework.h"
+
+Vec2 applyHomography(const Homography& H, const Vec2& m) {
+ double d = H[6] * m[0] + H[7] * m[1] + 1;
+ return Vec2((H[0] * m[0] + H[1] * m[1] + H[2]) / d,
+ (H[3] * m[0] + H[4] * m[1] + H[5]) / d);
+}
+
+Homography homographyFrom4Pairs(const Vec2 m1[4], const Vec2 m2[4]) {
+ Matrix<double> A(8,8);
+ Vector<double> B(8);
+ A.fill(0);
+ for(int k = 0; k < 4; k++) {
+ A(2*k, 0) = m1[k][0];
+ A(2*k, 1) = m1[k][1];
+ A(2*k, 2) = 1;
+ A(2*k, 6) = - m2[k][0] * m1[k][0];
+ A(2*k, 7) = - m2[k][0] * m1[k][1];
+ A(2*k + 1, 3) = m1[k][0];
+ A(2*k + 1, 4) = m1[k][1];
+ A(2*k + 1, 5) = 1;
+ A(2*k + 1, 6) = - m2[k][1] * m1[k][0];
+ A(2*k + 1, 7) = - m2[k][1] * m1[k][1];
+ B[2*k ] = m2[k][0];
+ B[2*k + 1] = m2[k][1];
+ }
+ Matrix<double> C=inverse(A);
+ if(norm(C)==0) // non inversible
+ return Homography(0.);
+ Vector<double> H=C*B;
+ return Homography(H.data()); // Painlesss Vector -> FVector conversion
+}
+
+Homography meanHomography(const Array<Feat>& feats1,const Array<Feat>
+& feats2,const MatchList& matches) {
+ int nb=int(matches.size());
+ Matrix<double> A(2*nb,8);
+ Vector<double> B(2*nb);
+ A.fill(0);
+ int k=0;
+ for (MatchList::const_iterator it=matches.begin();it!=matches.end();k++,it++){
+ Vec2 m1=feats1[it->first].pos;
+ Vec2 m2=feats2[it->second].pos;
+ A(2*k, 0) = m1[0];
+ A(2*k, 1) = m1[1];
+ A(2*k, 2) = 1;
+ A(2*k, 6) = - m2[0] * m1[0];
+ A(2*k, 7) = - m2[0] * m1[1];
+ A(2*k + 1, 3) = m1[0];
+ A(2*k + 1, 4) = m1[1];
+ A(2*k + 1, 5) = 1;
+ A(2*k + 1, 6) = - m2[1] * m1[0];
+ A(2*k + 1, 7) = - m2[1] * m1[1];
+ B[2*k ] = m2[0];
+ B[2*k + 1] = m2[1];
+ }
+ Matrix<double> C=pinverse(A); // Moindres carrés
+ if(norm(C)==0) // non invertible
+ return Homography(0.);
+ Vector<double> H=C*B;
+ return Homography(H.data()); // Painlesss Vector -> FVector conversion
+}
+
+Homography manualHomography(int w) {
+ Vec2 m1[4],m2[4];
+ int x, y;
+ for (int i = 0; i < 4; i++) {
+ GetMouse(x, y);
+ DrawCircle(x, y, 2, Color(255, 0, 0), 2);
+ if (x < w)
+ m1[i] = Vec2(x, y);
+ else
+ m2[i] = Vec2(x - w, y);
+ GetMouse(x, y);
+ DrawCircle(x, y, 2, Color(255, 0, 0), 2);
+ if (x < w)
+ m1[i] = Vec2(x, y);
+ else
+ m2[i] = Vec2(x - w, y);
+ }
+ return homographyFrom4Pairs(m1,m2);
+}
+
+class HomEstimator {
+ const Array<Feat>& feats1,feats2;
+public:
+ HomEstimator(const Array<Feat>& ft1,const Array<Feat>& ft2):feats1(ft1),feats2(ft2){}
+ template <class H_IT>
+ void operator() (const FArray<MatchList::iterator,4>& pairs, H_IT homographies) const {
+ Vec2 m1[4],m2[4];
+ for(int k=0;k<4;k++){
+ m1[k]=feats1[pairs[k]->first].pos;
+ m2[k]=feats2[pairs[k]->second].pos;
+ }
+ Homography H = homographyFrom4Pairs(m1, m2);
+ if (norm(H)!=0)
+ *homographies++ = H;
+ }
+};
+
+// Functor: estimates the square residual of a pair w.r.t. an Homography
+class HomResidual {
+ const Array<Feat>& feats1,feats2;
+public:
+ HomResidual(const Array<Feat>& ft1,const Array<Feat>& ft2):feats1(ft1),feats2(ft2){}
+ double operator() (const Homography& H,const pair<size_t,size_t>& p) const {
+ Vec2 m1(feats1[p.first].pos);
+ Vec2 m2(feats2[p.second].pos);
+ return norm2(applyHomography(H, m1) - m2);
+ }
+};
+
+void Ransac(const Array<Feat> & feats1, const Array<Feat> & feats2,
+ MatchList & matches) {
+ Homography H;
+ double outlier_thres;
+ double median_res = LMedS<4>(matches.begin(), matches.end(),
+ HomEstimator(feats1,feats2),
+ HomResidual(feats1,feats2),
+ H, &outlier_thres);
+ MatchList inliers;
+ for (MatchList::const_iterator it = matches.begin(); it != matches.end(); it++) {
+ if (HomResidual(feats1, feats2)(H, *it) < outlier_thres)
+ inliers.push_front(*it);
+ }
+ matches = inliers;
+}
diff --git a/framework.h b/framework.h
new file mode 100644
index 0000000..230a594
--- /dev/null
+++ b/framework.h
@@ -0,0 +1,32 @@
+#ifndef __FRAMEWORK_H_191109__
+#define __FRAMEWORK_H_191109__
+
+#include <iostream>
+#include <Imagine/Images.h>
+#include <Imagine/LinAlg.h>
+
+#include "Imagine/Features.h"
+#include "Imagine/Optim.h"
+
+using namespace std;
+using namespace Imagine;
+
+typedef FVector<double, 2> Vec2;
+typedef FVector<double, 3> Vec3;
+typedef FVector<double, 8> Homography;
+//TODO un truc moins immonde
+#ifndef FEAT_DE
+#define FEAT_DE
+typedef SIFTDetector Detector;
+typedef Detector::Feature Feat;
+#endif
+
+Vec2 applyHomography(const Homography& H, const Vec2& m);
+Homography homographyFrom4Pairs(const Vec2 m1[4], const Vec2 m2[4]);
+Homography meanHomography(const Array<Feat>& feats1,const Array<Feat>
+ & feats2,const MatchList& matches);
+Homography manualHomography(int w);
+void Ransac (const Array<Feat> & feats1, const Array<Feat> & feats2,
+ MatchList & matches);
+
+#endif