---
milena/sandbox/anthony/Makefile | 8 +-
milena/sandbox/anthony/new-implem/Makefile | 24 -
milena/sandbox/anthony/new-implem/src/keypoint.hh | 48 --
milena/sandbox/anthony/new-implem/src/main.cc | 72 --
milena/sandbox/anthony/new-implem/src/matrix.cc | 154 ----
milena/sandbox/anthony/new-implem/src/matrix.hh | 37 -
.../anthony/{new-implem => }/src/dog-space.hh | 0
milena/sandbox/anthony/src/keypoint.cc | 45 --
milena/sandbox/anthony/src/keypoint.hh | 58 +-
milena/sandbox/anthony/src/main.cc | 817 ++------------------
milena/sandbox/anthony/src/matrix.hh | 3 +-
.../anthony/{new-implem => }/src/scale-space.hh | 0
.../sandbox/anthony/{new-implem => }/src/vrac.hh | 0
13 files changed, 86 insertions(+), 1180 deletions(-)
delete mode 100644 milena/sandbox/anthony/new-implem/Makefile
delete mode 100644 milena/sandbox/anthony/new-implem/src/keypoint.hh
delete mode 100644 milena/sandbox/anthony/new-implem/src/main.cc
delete mode 100644 milena/sandbox/anthony/new-implem/src/matrix.cc
delete mode 100644 milena/sandbox/anthony/new-implem/src/matrix.hh
rename milena/sandbox/anthony/{new-implem => }/src/dog-space.hh (100%)
delete mode 100644 milena/sandbox/anthony/src/keypoint.cc
rename milena/sandbox/anthony/{new-implem => }/src/scale-space.hh (100%)
rename milena/sandbox/anthony/{new-implem => }/src/vrac.hh (100%)
diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/Makefile
index 43d90fd..5441d9a 100644
--- a/milena/sandbox/anthony/Makefile
+++ b/milena/sandbox/anthony/Makefile
@@ -4,12 +4,12 @@ CFLAGS=-Wall -Wextra
CLIBS=-I../../
CLEAN=*.o output/* log debug_*
-SRC=src/main.cc src/matrix.cc src/keypoint.cc
+SRC=src/main.cc src/matrix.cc
OUTPUT=a.out
-all: scale
+all: sift
-scale: clean
+sift: clean
$(CCACHE) $(CC) $(CFLAGS) -O3 -DNDEBUG $(CLIBS) $(SRC) -o $(OUTPUT)
debug: mrproper
@@ -21,4 +21,4 @@ clean:
mrproper: clean
rm -f $(OUTPUT)
-.PHONY: scale clean mrproper
+.PHONY: sift clean mrproper
diff --git a/milena/sandbox/anthony/new-implem/Makefile b/milena/sandbox/anthony/new-implem/Makefile
deleted file mode 100644
index a4f63d9..0000000
--- a/milena/sandbox/anthony/new-implem/Makefile
+++ /dev/null
@@ -1,24 +0,0 @@
-CCACHE=
-CC=g++
-CFLAGS=-Wall -Wextra
-CLIBS=-I../../../
-CLEAN=*.o output/* log debug_*
-
-SRC=src/main.cc src/matrix.cc
-OUTPUT=new-implem
-
-all: new-implem
-
-new-implem: clean
- $(CCACHE) $(CC) $(CFLAGS) -O3 -DNDEBUG $(CLIBS) $(SRC) -o $(OUTPUT)
-
-debug: mrproper
- $(CCACHE) $(CC) $(CFLAGS) -g -ggdb $(CLIBS) $(SRC) -o $(OUTPUT)
-
-clean:
- rm -rf $(CLEAN)
-
-mrproper: clean
- rm -f $(OUTPUT)
-
-.PHONY: new-implem clean mrproper
diff --git a/milena/sandbox/anthony/new-implem/src/keypoint.hh b/milena/sandbox/anthony/new-implem/src/keypoint.hh
deleted file mode 100644
index b6b2415..0000000
--- a/milena/sandbox/anthony/new-implem/src/keypoint.hh
+++ /dev/null
@@ -1,48 +0,0 @@
-#ifndef KEYPOINT_HH
-# define KEYPOINT_HH
-
-# include <ostream>
-
-class Keypoint
-{
- public:
- Keypoint(int _x, int _y,
- unsigned _octave, unsigned _gradient,
- bool _maximum)
- {
- x = _x;
- y = _y;
- octave = _octave;
- gradient = _gradient;
- maximum = _maximum;
- }
-
- inline int getX() const { return x; }
- inline int getY() const { return y; }
- inline unsigned getOctave() const { return octave; }
- inline unsigned getGradient() const { return gradient; }
- inline bool getType() const { return maximum; }
- inline float getOrientation() const { return orientation; }
-
- inline void setX(int _x) { x = _x; }
- inline void setY(int _y) { y = _y; }
- inline void setOrientation(float _o) { orientation = _o; }
-
- friend std::ostream& operator<< (std::ostream& stream, const Keypoint& k)
- {
- stream << "(" << k.getX() << "," << k.getY() << ")"
- << " o:" << k.getOctave() << " g:" << k.getGradient();
-
- return stream;
- }
-
- private:
- int x;
- int y;
- unsigned octave;
- unsigned gradient;
- bool maximum;
- float orientation;
-};
-
-#endif /* ! KEYPOINT_HH */
diff --git a/milena/sandbox/anthony/new-implem/src/main.cc b/milena/sandbox/anthony/new-implem/src/main.cc
deleted file mode 100644
index afbb476..0000000
--- a/milena/sandbox/anthony/new-implem/src/main.cc
+++ /dev/null
@@ -1,72 +0,0 @@
-#include <mln/core/image/image2d.hh>
-#include <mln/data/all.hh>
-#include <mln/io/pgm/all.hh>
-#include <mln/io/ppm/save.hh>
-#include <mln/core/alias/neighb2d.hh>
-#include <mln/literal/all.hh>
-#include <mln/value/int_u8.hh>
-#include <mln/value/int_s16.hh>
-
-//#include <cmath>
-#include <string>
-#include <list>
-
-#include "scale-space.hh"
-#include "dog-space.hh"
-#include "keypoint.hh"
-//#include "matrix.hh"
-
-int main(void)
-{
- typedef image2d<value::int_u8> I;
- typedef image2d<value::int_s16> S;
- typedef image2d<value::rgb8> C;
-
- // General parameters
- std::string source("images/keith.pbm");
- const unsigned blur_level = 5;
- const unsigned octave_level = 4;
-
- ScaleSpace<I> *ss = new ScaleSpace<I>(octave_level, blur_level);
- DoGSpace<I, S>* dogs = new DoGSpace<I, S>(ss);
- std::list<Keypoint> keypoints;
-
- I original;
- C extrema, extrema_b,
- extrema2, extrema2_b;
- //improved, improved_b;
-
- io::pgm::load(original, source.c_str());
-
- initialize(extrema_b, original);
- initialize(extrema2_b, original);
- //initialize(improved_b, original);
-
- extrema = data::convert(value::rgb8(), original);
- extrema2 = data::convert(value::rgb8(), original);
- //improved = data::convert(value::rgb8(), original);
-
- // Localization
- ss->build(original);
- //ss->save();
- dogs->build();
- //dogs->save();
- dogs->findKeypoints(keypoints);
-
- writeKeypoints(keypoints, extrema, extrema_b);
-
- dogs->discardLowContrastKeypoints(keypoints);
- dogs->eliminateEdgeResponses(keypoints);
-
- writeKeypoints(keypoints, extrema2, extrema2_b);
-
- // Save
- io::ppm::save(extrema, "output/extrema.ppm");
- io::ppm::save(extrema_b, "output/extrema_b.ppm");
- io::ppm::save(extrema2, "output/extrema2.ppm");
- io::ppm::save(extrema2_b, "output/extrema2_b.ppm");
- //io::ppm::save(improved, "output/extrema_improved.ppm");
- //io::ppm::save(improved_b, "output/extrema_b_improved.ppm");
-
- return 0;
-}
diff --git a/milena/sandbox/anthony/new-implem/src/matrix.cc b/milena/sandbox/anthony/new-implem/src/matrix.cc
deleted file mode 100644
index ee00ffc..0000000
--- a/milena/sandbox/anthony/new-implem/src/matrix.cc
+++ /dev/null
@@ -1,154 +0,0 @@
-#include "matrix.hh"
-
-Matrix::Matrix(unsigned h, unsigned w)
-{
- height = h;
- width = w;
-}
-
-void Matrix::initialize()
-{
- for (unsigned i = 0; i < height; ++i)
- {
- std::vector<float> line;
-
- for (unsigned j = 0; j < width; ++j)
- line.push_back(0);
-
- m.push_back(line);
- }
-}
-
-void Matrix::initialize(float *values)
-{
- for (unsigned i = 0; i < height; ++i)
- {
- std::vector<float> line;
-
- for (unsigned j = 0; j < width; ++j)
- line.push_back(values[i * width + j]);
-
- m.push_back(line);
- }
-}
-
-void Matrix::print()
-{
- for (unsigned i = 0; i < height; ++i)
- {
- std::cout << "|";
-
- for (unsigned j = 0; j < width; ++j)
- std::cout << " " << m.at(i).at(j);
-
- std::cout << " |" << std::endl;
- }
- std::cout << std::endl;
-}
-
-float Matrix::determinant()
-{
- float d1 = m.at(0).at(0) * m.at(1).at(1) * m.at(2).at(2);
- float d2 = m.at(0).at(1) * m.at(1).at(2) * m.at(2).at(0);
- float d3 = m.at(0).at(2) * m.at(1).at(0) * m.at(2).at(1);
-
- float d4 = m.at(0).at(2) * m.at(1).at(1) * m.at(2).at(0);
- float d5 = m.at(1).at(2) * m.at(2).at(1) * m.at(0).at(0);
- float d6 = m.at(2).at(2) * m.at(0).at(1) * m.at(1).at(0);
-
- return (d1 + d2 + d3 - d4 - d5 - d6);
-}
-
-bool Matrix::inverse(Matrix& inversed)
-{
- bool canInverse = false;
-
- float det = determinant();
-
- if (det != 0)
- {
- float a = m.at(0).at(0);
- float b = m.at(0).at(1);
- float c = m.at(0).at(2);
-
- float d = m.at(1).at(0);
- float e = m.at(1).at(1);
- float f = m.at(1).at(2);
-
- float g = m.at(2).at(0);
- float h = m.at(2).at(1);
- float i = m.at(2).at(2);
-
- Matrix comatrix(3, 3);
- comatrix.initialize();
-
- comatrix.set(0, 0, e*i - f*h);
- comatrix.set(0, 1, c*h - b*i);
- comatrix.set(0, 2, b*f - c*e);
-
- comatrix.set(1, 0, f*g - d*i);
- comatrix.set(1, 1, a*i - c*g);
- comatrix.set(1, 2, c*d - a*f);
-
- comatrix.set(2, 0, d*h - e*g);
- comatrix.set(2, 1, b*g - a*h);
- comatrix.set(2, 2, a*e - b*d);
-
- for (unsigned i = 0; i < comatrix.getWidth(); ++i)
- for (unsigned j = 0; j < comatrix.getHeight(); ++j)
- comatrix.set(i, j, comatrix.get(i,j) / det);
-
- inversed = comatrix;
-
- canInverse = true;
- }
-
- return canInverse;
-}
-
-// Operators
-Matrix Matrix::operator*(Matrix& right)
-{
- Matrix result(height, right.getWidth());
- result.initialize();
-
- for (unsigned i = 0; i < height; ++i)
- for (unsigned j = 0; j < width; ++j)
- for (unsigned k = 0; k < right.getWidth(); ++k)
- result.add(i, k, m.at(i).at(j) * right.get(j, k));
-
- return result;
-}
-
-void Matrix::absolute()
-{
- for (unsigned i = 0; i < height; ++i)
- for (unsigned j = 0; j < width; ++j)
- m.at(i).at(j) *= (m.at(i).at(j) >= 0) ? 1 : -1;
-}
-
-bool Matrix::isMajorOffset()
-{
- unsigned i = 0;
- unsigned j = 0;
- bool isMajorOffset = false;
-
- while (!isMajorOffset && i < height)
- {
- while (!isMajorOffset && j < width)
- {
- isMajorOffset = m.at(i).at(j) > 0.5 || - m.at(i).at(j) > 0.5;
- ++j;
- }
- ++i;
- }
-
- return isMajorOffset;
-}
-
-// Getters and Setters
-unsigned Matrix::getHeight() { return height; }
-unsigned Matrix::getWidth() { return width; }
-float Matrix::get(unsigned i, unsigned j) { return m.at(i).at(j); }
-void Matrix::set(unsigned i, unsigned j, float x) { m.at(i).at(j) = x; }
-void Matrix::add(unsigned i, unsigned j, float x) { m.at(i).at(j) += x; }
diff --git a/milena/sandbox/anthony/new-implem/src/matrix.hh b/milena/sandbox/anthony/new-implem/src/matrix.hh
deleted file mode 100644
index c42de59..0000000
--- a/milena/sandbox/anthony/new-implem/src/matrix.hh
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifndef MATRIX_HH
-# define MATRIX_HH
-
-# include <iostream>
-# include <vector>
-
-class Matrix
-{
- public:
- // Constructor
- Matrix(unsigned i, unsigned j);
-
- void initialize();
- void initialize(float *values);
- void print();
- bool inverse(Matrix& inversed);
- void absolute();
- bool isMajorOffset();
- float determinant();
-
- // Operators
- Matrix operator*(Matrix& right);
-
- // Getters and Setters
- unsigned getHeight();
- unsigned getWidth();
- float get(unsigned i, unsigned j);
- void set(unsigned i, unsigned j, float value);
- void add(unsigned i, unsigned j, float value);
-
- private:
- std::vector< std::vector<float> > m;
- unsigned height;
- unsigned width;
-};
-
-#endif /* ! MATRIX_HH */
diff --git a/milena/sandbox/anthony/new-implem/src/dog-space.hh b/milena/sandbox/anthony/src/dog-space.hh
similarity index 100%
rename from milena/sandbox/anthony/new-implem/src/dog-space.hh
rename to milena/sandbox/anthony/src/dog-space.hh
diff --git a/milena/sandbox/anthony/src/keypoint.cc b/milena/sandbox/anthony/src/keypoint.cc
deleted file mode 100644
index 937e78d..0000000
--- a/milena/sandbox/anthony/src/keypoint.cc
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "keypoint.hh"
-
-// Constructor
-Keypoint::Keypoint(int i, int j,
- int scale, int size, bool maximum)
-{
- this->i = i;
- this->j = j;
- this->scale = scale;
- this->size = size;
- this->maximum = maximum;
-}
-
-bool Keypoint::add(Matrix& offset, Keypoint& old)
-{
- if (offset.get(0, 0) >= 1 || offset.get(0, 0) <= -1
- || offset.get(1, 0) >= 1 || offset.get(1, 0) <= -1)
- {
- i += offset.get(0, 0);
- j += offset.get(1, 0);
- scale += offset.get(2, 0);
-
- return (i != old.getI() && j != old.getJ());
- }
- else
- return false;
-}
-
-//mln::point2d Keypoint::getPoint()
-//{
- //mln::point2d p(i, j);
-
- //return p;
-//}
-
-// Getters
-int Keypoint::getI() { return i; }
-int Keypoint::getJ() { return j; }
-int Keypoint::getScale() { return scale; }
-int Keypoint::getSize() { return size; }
-
-// Setters
-void Keypoint::setI(int offsetI) { i = offsetI; }
-void Keypoint::setJ(int offsetJ) { j = offsetJ; }
-void Keypoint::setScale(int offsetScale) { scale = offsetScale; }
diff --git a/milena/sandbox/anthony/src/keypoint.hh b/milena/sandbox/anthony/src/keypoint.hh
index af24a4c..b6b2415 100644
--- a/milena/sandbox/anthony/src/keypoint.hh
+++ b/milena/sandbox/anthony/src/keypoint.hh
@@ -1,36 +1,48 @@
#ifndef KEYPOINT_HH
# define KEYPOINT_HH
-//# include <mln/core/alias/point2d.hh>
# include <ostream>
-# include "matrix.hh"
class Keypoint
{
public:
- // Constructor
- Keypoint(int i, int j,
- int scale, int size, bool maximum);
-
- bool add(Matrix& k, Keypoint& old);
- //mln::point2d getPoint();
-
- // Getters
- int getI();
- int getJ();
- int getScale();
- int getSize();
-
- void setI(int offsetI);
- void setJ(int offsetJ);
- void setScale(int offsetScale);
+ Keypoint(int _x, int _y,
+ unsigned _octave, unsigned _gradient,
+ bool _maximum)
+ {
+ x = _x;
+ y = _y;
+ octave = _octave;
+ gradient = _gradient;
+ maximum = _maximum;
+ }
+
+ inline int getX() const { return x; }
+ inline int getY() const { return y; }
+ inline unsigned getOctave() const { return octave; }
+ inline unsigned getGradient() const { return gradient; }
+ inline bool getType() const { return maximum; }
+ inline float getOrientation() const { return orientation; }
+
+ inline void setX(int _x) { x = _x; }
+ inline void setY(int _y) { y = _y; }
+ inline void setOrientation(float _o) { orientation = _o; }
+
+ friend std::ostream& operator<< (std::ostream& stream, const Keypoint& k)
+ {
+ stream << "(" << k.getX() << "," << k.getY() << ")"
+ << " o:" << k.getOctave() << " g:" << k.getGradient();
+
+ return stream;
+ }
private:
- int i;
- int j;
- int scale;
- int size;
- bool maximum;
+ int x;
+ int y;
+ unsigned octave;
+ unsigned gradient;
+ bool maximum;
+ float orientation;
};
#endif /* ! KEYPOINT_HH */
diff --git a/milena/sandbox/anthony/src/main.cc b/milena/sandbox/anthony/src/main.cc
index 59fa1ce..afbb476 100644
--- a/milena/sandbox/anthony/src/main.cc
+++ b/milena/sandbox/anthony/src/main.cc
@@ -4,794 +4,69 @@
#include <mln/io/ppm/save.hh>
#include <mln/core/alias/neighb2d.hh>
#include <mln/literal/all.hh>
-#include <mln/value/int_s8.hh>
+#include <mln/value/int_u8.hh>
+#include <mln/value/int_s16.hh>
-#include <cmath>
+//#include <cmath>
+#include <string>
+#include <list>
+#include "scale-space.hh"
+#include "dog-space.hh"
#include "keypoint.hh"
-#include "matrix.hh"
+//#include "matrix.hh"
-using namespace mln;
-
-// Set the gaussian kernel
-float *gaussian_kernel()
-{
- float variance = sqrt(2.0);
- unsigned index = 0;
- float *kernel = new float[49];
- static const float pi = 4.0 * atan(1.0);
- float div = (2.0 * pi * variance);
-
- for (int i = -3; i < 4; ++i)
- {
- for (int j = -3; j < 4; ++j)
- {
- float e = exp(- (i*i + j*j) / (2.0 * variance) );
- kernel[index++] = e / div;
- }
- }
-
- return kernel;
-}
-
-// Apply the convolution of the image by the kernel
-template<typename I>
-void convolve(const I& original, I& filtered, float *kernel)
-{
- mln_piter(I) p(filtered.domain());
- window2d win;
-
- for (int i = -3; i < 4; ++i)
- for (int j = -3; j < 4; ++j)
- win.insert(i, j);
-
- // Iterate through all image sites
- for_all(p)
- {
- // Create the window around the site
- mln_qiter(window2d) q(win, p);
- float sum = 0;
- int index = 0;
-
- // Iterate through all window image sites
- for_all(q)
- {
- if (filtered.has(q))
- sum += original(q) * kernel[index++];
- else
- sum += 127 * kernel[index++];
- }
-
- filtered(p) = static_cast<int>(sum);
- }
-}
-
-// Blur the image thanks to a convolution matrix
-template<typename I>
-I blur(const I& original)
-{
- I filtered;
- float *kernel = gaussian_kernel();
-
- initialize(filtered, original);
- convolve(original, filtered, kernel);
-
- return filtered;
-}
-
-// Downscale by 2 the image resolution
-template<typename I>
-I downscaling2(const I& original)
-{
- mln_piter(I) p(original.domain());
- window2d win;
- int index = 0;
-
- // Initialize the rescaled image
- box2d size(original.nrows() / 2, original.ncols() / 2);
- I reduced(size);
-
- // Set the 2x2 window matrix
- for (int i = 0; i < 2; ++i)
- for (int j = 0; j < 2; ++j)
- win.insert(i, j);
-
- // Iterate through all image sites
- for_all(p)
- {
- // Get 1/4 of total pixels
- if (index % 2 == 0)
- {
- int i = (index % original.ncols()) / 2;
- int j = (index / original.ncols()) / 2;
- int count = 0;
- int sum = 0;
- mln_qiter(window2d) q(win, p);
-
- // Iterate through all window image sites
- for_all(q)
- {
- if (original.has(q))
- {
- sum += original(q);
- ++count;
- }
- }
-
- if (reduced.has(point2d(j, i)))
- opt::at(reduced, j, i) = sum / count;
- }
-
- ++index;
- }
-
- return reduced;
-}
-
-float *upscalingTL_kernel()
-{
- const static float one = 1.0 / 8.0;
- const static float two = 2.0 / 8.0;
- const static float three = 2.0 / 8.0;
- const static float four = 3.0 / 8.0;
-
- float *kernel = new float[4];
-
- kernel[0] = one;
- kernel[1] = two;
- kernel[2] = three;
- kernel[3] = four;
-
- return kernel;
-}
-
-float *upscalingTR_kernel()
-{
- const static float one = 2.0 / 8.0;
- const static float two = 1.0 / 8.0;
- const static float three = 3.0 / 8.0;
- const static float four = 2.0 / 8.0;
-
- float *kernel = new float[4];
-
- kernel[0] = one;
- kernel[1] = two;
- kernel[2] = three;
- kernel[3] = four;
-
- return kernel;
-}
-
-float *upscalingBL_kernel()
-{
- const static float one = 2.0 / 8.0;
- const static float two = 3.0 / 8.0;
- const static float three = 1.0 / 8.0;
- const static float four = 2.0 / 8.0;
-
- float *kernel = new float[4];
-
- kernel[0] = one;
- kernel[1] = two;
- kernel[2] = three;
- kernel[3] = four;
-
- return kernel;
-}
-
-float *upscalingBR_kernel()
-{
- const static float one = 3.0 / 8.0;
- const static float two = 2.0 / 8.0;
- const static float three = 2.0 / 8.0;
- const static float four = 1.0 / 8.0;
-
- float *kernel = new float[4];
-
- kernel[0] = one;
- kernel[1] = two;
- kernel[2] = three;
- kernel[3] = four;
-
- return kernel;
-}
-
-template<typename I>
-unsigned upscalePixel(const I& original, int i, int j, int type)
-{
- float *kernel;
- window2d win;
- mln_psite(I) p(i, j);
- float value = 0;
- unsigned index = 0;
- int iStart = 0;
- int jStart = 0;
-
- if (type == 0)
- {
- kernel = upscalingTL_kernel();
- iStart = -1;
- jStart = -1;
- }
- else if (type == 1)
- {
- kernel = upscalingTR_kernel();
- iStart = 0;
- jStart = -1;
- }
- else if (type == 2)
- {
- kernel = upscalingBL_kernel();
- iStart = -1;
- jStart = 0;
- }
- else
- {
- kernel = upscalingBR_kernel();
- iStart = 0;
- jStart = 0;
- }
-
- for (int ii = iStart; ii < iStart + 2; ++ii)
- for (int jj = jStart; jj < jStart + 2; ++jj)
- win.insert(ii, jj);
-
- mln_qiter(window2d) q(win, p);
-
- for_all(q)
- {
- if (original.has(q))
- {
- value += original(q) * kernel[index];
- }
-
- ++index;
- }
-
- return (static_cast<int>(value));
-}
-
-// Upscale by 2 the image resolution
-template<typename I>
-I upscaling2(const I& original)
-{
- box2d size(original.nrows() * 2, original.ncols() * 2);
- I upscaled(size);
-
- mln_piter(I) p(upscaled.domain());
- int index = 0;
-
- for_all(p)
- {
- int i = p.col();
- int j = p.row();
-
- if (i % 2 == 0 && j % 2 == 0) // Top left
- opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 0);
- else if (i % 2 == 0) // Bottom left
- opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 2);
- else if (j % 2 == 0) // Top right
- opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 1);
- else // Bottom right
- opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 3);
-
- ++index;
- }
-
-
- return upscaled;
-}
-
-// Build the n-th octave on b blur levels
-template<typename I>
-void buildOctave(const I& original,
- unsigned n,
- unsigned b,
- std::vector< std::vector<I> >& scaleSpace)
-{
- std::vector<I> octave;
- std::stringstream name;
- I blured;
- initialize(blured, original);
-
- name << "output/SS-o" << n << "-0.pgm";
- io::pgm::save(original, name.str());
-
- octave.push_back(original);
- blured = duplicate(original);
-
- for(unsigned i = 1; i <= b; i *= 2)
- {
- blured = blur(blured);
- octave.push_back(blured);
-
- name.str("");
- name << "output/SS-o" << n << "-" << i << ".pgm";
- io::pgm::save(blured, name.str());
- }
-
- scaleSpace.push_back(octave);
-}
-
-// Build n octaves on b blur levels
-template<typename I>
-void buildScaleSpace(std::vector< std::vector<I> >& scaleSpace,
- const I& original,
- unsigned n,
- unsigned b)
-{
- if (n > 2)
- {
- // Upscaled octave
- I upscaled = upscaling2(original);
- buildOctave(upscaled, 0, b, scaleSpace);
-
- // 1st octave
- buildOctave(original, 1, b, scaleSpace);
-
- // 2nd octave
- I downscaled = downscaling2(original);
- buildOctave(downscaled, 2, b, scaleSpace);
-
- for (unsigned i = 3; i <= n; ++i)
- {
- // i-th octave
- downscaled = downscaling2(downscaled);
- buildOctave(downscaled, i, b, scaleSpace);
- }
- }
- else
- {
- std::cerr << "ERROR: buildScaleSpace: Wrong level number (" << n << ")"
- << std::endl;
- }
-}
-
-// Build the DoG pyramid of the scale space
-template<typename I>
-void buildDifferenceOfGaussianSpace(std::vector< std::vector<I> >& scaleSpace,
- std::vector< std::vector<I> >& dogSpace)
-{
- for (unsigned i = 0; i < scaleSpace.size(); ++i)
- {
- std::vector<I> octave = scaleSpace.at(i);
- std::vector<I> dogOctave;
-
- for (unsigned j = 0; j < octave.size() - 1; ++j)
- {
- std::stringstream name;
- I level1 = octave.at(j);
- I level2 = octave.at(j + 1);
- I dog;
-
- initialize(dog, level1);
-
- mln_piter(I) p(level1.domain());
-
- for_all(p)
- {
- int diff = ((int) level1(p)) - ((int) level2(p));
- dog(p) = (diff >= 0) ? diff : 0;
- }
-
- dogOctave.push_back(dog);
- name << "output/DoG_o" << i << "-" << j << ".pgm";
- io::pgm::save(dog, name.str());
- }
-
- dogSpace.push_back(dogOctave);
- }
-}
-
-// Find the extremum value in the c8 neighborhood
-template<typename T, typename I, typename N>
-void find_extremum(T& min, T& max, const I& current, N n)
-{
- for_all(n)
- {
- if (current(n) > max)
- max = current(n);
-
- if (current(n) < min)
- min = current(n);
- }
-}
-
-// Find the minimum value in the upper and lower layers around p
-template<typename T, typename I, typename P>
-void find_min(T& min, const I& upper, const I& lower, const P& pUpper, const P& pLower)
-{
- mln_niter(neighb2d) nUpper(c8(), pUpper);
- mln_niter(neighb2d) nLower(c8(), pLower);
- T tmp = min;
-
- // Upper processing
- for_all(nUpper)
- if (upper(nUpper) < min)
- min = upper(nUpper);
-
- // Lower processing
- for_all(nLower)
- if (lower(nLower) < min)
- min = lower(nLower);
-
- std::cout << tmp << " >= " << min << std::endl;
-}
-
-// Find the maximum value in the upper and lower layers around p
-template<typename T, typename I, typename P>
-void find_max(T& max, const I& upper, const I& lower, const P& pUpper, const P& pLower)
-{
- mln_niter(neighb2d) nUpper(c8(), pUpper);
- mln_niter(neighb2d) nLower(c8(), pLower);
- T tmp = max;
-
- // Upper processing
- for_all(nUpper)
- if (upper(nUpper) > max)
- max = upper(nUpper);
-
- // Lower processing
- for_all(nLower)
- if (lower(nLower) > max)
- max = lower(nLower);
-
- std::cout << tmp << " >= " << max << std::endl;
-}
-
-template<typename T, typename I, typename P>
-bool hasEdgeResponse(const I& dog, const P& p, const T& center)
-{
- const static float bound_ratio = 12.1;
-
- /* Set the neighbors:
- * n5 n1 n6
- * n2 p n3
- * n7 n4 n8
- */
- float n1 = dog(mln_psite(I)(p.row() - 1, p.col()));
- float n2 = dog(mln_psite(I)(p.row(), p.col() - 1));
- float n3 = dog(mln_psite(I)(p.row(), p.col() + 1));
- float n4 = dog(mln_psite(I)(p.row() + 1, p.col()));
-
- float n5 = dog(mln_psite(I)(p.row() - 1, p.col() - 1));
- float n6 = dog(mln_psite(I)(p.row() - 1, p.col() + 1));
- float n7 = dog(mln_psite(I)(p.row() + 1, p.col() - 1));
- float n8 = dog(mln_psite(I)(p.row() + 1, p.col() + 1));
-
- // Set Hessian matrix and find the ratio
-
- float dxx = n3 - 2.0 * center + n2;
- float dyy = n4 - 2.0 * center + n1;
-
- float dxy = (n8 - n6 - n7 + n5) / 4.0;
- float trace = dxx + dyy;
- float determinant = dxx * dyy - dxy * dxy;
-
- if (determinant <= 0)
- return true;
-
- float ratio = (trace * trace) / determinant;
-
- return (ratio >= bound_ratio);
-}
-
-// Build the extrema image
-template<typename I, typename C>
-void buildExtrema(C& extrema,
- std::vector< std::vector<I> >& dogSpace,
- std::vector<Keypoint>& keypoints)
-{
- // Number of scales (gaussian meaning)
- unsigned scales = dogSpace.at(0).size();
- std::cout << "Scales: " << scales << std::endl
- << "Resolutions: " << dogSpace.size() << std::endl;
-
-
- // Iterates through all scales
- for (unsigned i = 0; i < scales - 2; ++i)
- {
- std::cout << "scales ENTER" << std::endl;
-
- // Iterates through all resolutions
- for (unsigned j = 1; j < dogSpace.size() - 2; ++j)
- {
- std::cout << "resolutions ENTER" << std::endl;
- I current = dogSpace.at(j).at(i);
- I upper = dogSpace.at(j+1).at(i);
- I lower = dogSpace.at(j-1).at(i);
-
- C debugDraw1 = data::convert(value::rgb8(), upper);
- C debugDraw2 = data::convert(value::rgb8(), lower);
-
-
- mln_piter(I) p(current.domain());
- int first = 0;
- int bound = 1000;
- unsigned threshold = 0;
-
-
- for_all(p)
- {
- value::int_u8 min = 255;
- value::int_u8 max = 0;
- value::int_u8 center = 0;
- ++first;
-
- mln_niter(neighb2d) n(c8(), p);
- center = current(p);
-
- mln_psite(I) pUpper(p.row() / 2, p.col() / 2);
- mln_psite(I) pLower(p.row() * 2, p.col() * 2);
-
- // Find the min and max in the current neighborhood (c8)
- find_extremum(min, max, current, n);
-
- if ((unsigned) center + threshold < min)
- {
-
- // Find the min in other levels of the pyramid
- find_min(min, upper, lower, pUpper, pLower);
-
- if ((unsigned) center + threshold < min && !hasEdgeResponse(current, p, center))
- {
- std::cout << center << " < " << min << std::endl;
-
- if (first > bound)
- {
- debugDraw1(pUpper) = literal::yellow;
- debugDraw2(pLower) = literal::yellow;
-
-
- ++first;
- }
- mln_psite(I) pOriginal(p.row() * pow(2, (j-1)),
- p.col() * pow(2, (j-1)));
- extrema(pOriginal) = literal::red;
- keypoints.push_back(Keypoint(p.row(), p.col(), i, j, false));
- }
- }
- else if (center > (unsigned) max + threshold)
- {
- // Find the max in other levels of the pyramid
- find_max(max, upper, lower, pUpper, pLower);
-
- if (center > (unsigned) max + threshold && !hasEdgeResponse(current, p, center))
- {
- std::cout << center << " > " << max << std::endl;
-
- if (first > bound)
- {
- debugDraw1(pUpper) = literal::yellow;
- debugDraw2(pLower) = literal::yellow;
-
- ++first;
- }
-
- mln_psite(I) pOriginal(p.row() * pow(2, (j-1)),
- p.col() * pow(2, (j-1)));
- extrema(pOriginal) = literal::red;
- keypoints.push_back(Keypoint(p.row(), p.col(), i, j, true));
- }
- }
- }
- io::ppm::save(debugDraw1, "debug_upper.ppm");
- io::ppm::save(debugDraw2, "debug_lower.ppm");
- }
- }
-}
-
-// Compute the first order matrix for interpolation
-template<typename I>
-void computeFirstOrderMatrix(const std::vector< std::vector<I> >& dogSpace,
- Keypoint& k,
- Matrix& matrix)
-{
- matrix.initialize();
-
- std::vector<I> current = dogSpace.at(1);
- std::vector<I> upper = dogSpace.at(2);
- std::vector<I> lower = dogSpace.at(0);
-
- /* Point naming convention
- * n1
- * n2 c n3
- * n4
- */
- point2d c(k.getI(), k.getJ());
-
- point2d n1(k.getI(), k.getJ() - 1);
- point2d n2(k.getI() - 1, k.getJ());
- point2d n3(k.getI() + 1, k.getJ());
- point2d n4(k.getI(), k.getJ() + 1);
-
- point2d cupper(c.row() / 2, c.col() / 2);
- point2d clower(c.row() * 2, c.col() * 2);
-
- // Compute partial x
- matrix.set(0, 0, (current.at(0)(n3) - current.at(0)(n2)) / 2.0);
-
- // Compute partial y
- matrix.set(1, 0, (current.at(0)(n3) - current.at(0)(n2)) / 2.0);
-
- // Compute partial scale
- matrix.set(2, 0, (lower.at(0)(clower) - upper.at(0)(cupper)) / 2.0);
-}
-
-// Compute the inversed second order matrix for interpolation
-template<typename I>
-void computeSecondOrderMatrix(const std::vector< std::vector<I> >& dogSpace,
- Keypoint& k,
- Matrix& matrix)
-{
- matrix.initialize();
-
- const static std::vector<I> current = dogSpace.at(1);
- const static I upper = dogSpace.at(2).at(0);
- const static I lower = dogSpace.at(0).at(0);
- const static I main = current.at(0);
-
- /* Point naming convention
- * n5 n1 n6
- * n2 c n3
- * n7 n4 n8
- */
- point2d c(k.getI(), k.getJ());
-
- point2d n1(k.getI(), k.getJ() - 1);
- point2d n2(k.getI() - 1, k.getJ());
- point2d n3(k.getI() + 1, k.getJ());
- point2d n4(k.getI(), k.getJ() + 1);
-
- point2d n5(k.getI() - 1, k.getJ() - 1);
- point2d n6(k.getI() + 1, k.getJ() - 1);
- point2d n7(k.getI() + 1, k.getJ() + 1);
- point2d n8(k.getI() - 1, k.getJ() + 1);
-
- point2d cupper(c.row() / 2, c.col() / 2);
- point2d n1upper(n1.row() / 2, n1.col() / 2);
- point2d n2upper(n2.row() / 2, n3.col() / 2);
- point2d n3upper(n3.row() / 2, n3.col() / 2);
- point2d n4upper(n4.row() / 2, n4.col() / 2);
-
- point2d clower(c.row() * 2, c.col() * 2);
- point2d n1lower(n1.row() * 2, n1.col() * 2);
- point2d n2lower(n2.row() * 2, n2.col() * 2);
- point2d n3lower(n3.row() * 2, n3.col() * 2);
- point2d n4lower(n4.row() * 2, n4.col() * 2);
-
- float partial_xx, partial_yy, partial_ss, partial_xy, partial_xs, partial_ys;
-
- // Compute partial xx
- partial_xx = main(n3) - 2 * main(c) + main(n2);
-
- // Compute partial yy
- partial_yy = main(n4) - 2 * main(c) + main(n1);
-
- // Compute partial ss
- partial_ss = lower(clower) - 2 * main(c)
- + upper(cupper);
-
- // Compute partial xy
- partial_xy = (main(n8) - main(n6) - main(n7) + main(n5)) / 4.0;
-
- // Compute partial xs
- partial_xs = (upper(n3upper) - lower(n3lower)
- - upper(n2upper) + lower(n2lower)) / 4.0;
-
- // Compute partial ys
- partial_ys = (upper(n4upper) - lower(n4lower)
- - upper(n1upper) + lower(n1lower)) / 4.0;
-
- matrix.set(0, 0, partial_xx);
- matrix.set(0, 1, partial_xy);
- matrix.set(0, 2, partial_xs);
-
- matrix.set(1, 0, partial_xy);
- matrix.set(1, 1, partial_yy);
- matrix.set(1, 2, partial_ys);
-
- matrix.set(2, 0, partial_xs);
- matrix.set(2, 1, partial_ys);
- matrix.set(2, 2, partial_ss);
-}
-
-// Discard low contrast keypoints thanks to taylor interpolation
-// See Brown and Lowe paper (2002)
-template<typename I, typename C>
-void discardLowContrastKeypoints(const std::vector< std::vector<I> >& dogSpace,
- std::vector<Keypoint>& keypoints,
- C& extrema)
-{
- Keypoint k_old(0, 0, 0, 0, true);
- bool hasMoved = false;
-
- for (unsigned i = 0; i < keypoints.size(); ++i)
- {
- Keypoint k = keypoints.at(i);
-
- Matrix firstOrder(3, 1);
- Matrix secondOrder(3, 3);
- Matrix inversed(secondOrder.getHeight(), secondOrder.getWidth());
-
- computeFirstOrderMatrix(dogSpace, k, firstOrder);
- computeSecondOrderMatrix(dogSpace, k, secondOrder);
- bool hasInversed = secondOrder.inverse(inversed);
-
- if (hasInversed)
- {
- Matrix x = inversed * firstOrder;
- bool isMajorOffset = x.isMajorOffset();
-
- if (isMajorOffset)
- {
- hasMoved = k.add(x, k_old);
-
- point2d p(k.getI(), k.getJ());
-
- if (hasMoved && extrema.has(p)
- && p.row() >= 0 && p.row() < geom::max_row(extrema)
- && p.col() >= 0 && p.col() < geom::max_col(extrema))
- {
- // extrema(p) = literal::red;
- k_old = keypoints.at(i);
- keypoints.at(i--) = k;
- }
- }
- else
- {
- point2d p(k.getI(), k.getJ());
-
- if (hasMoved)
- extrema(p) = literal::green;
- else
- extrema(p) = literal::red;
-
- hasMoved = false;
- }
- }
- }
-}
-
-// Main entry point
-int main(void)//int argc, char** argv)
+int main(void)
{
typedef image2d<value::int_u8> I;
+ typedef image2d<value::int_s16> S;
typedef image2d<value::rgb8> C;
+ // General parameters
+ std::string source("images/keith.pbm");
const unsigned blur_level = 5;
- const unsigned octave_level = 3;
- std::vector< std::vector<I> > scaleSpace;
- std::vector< std::vector<I> > dogSpace;
- std::vector<Keypoint> keypoints;
+ const unsigned octave_level = 4;
+
+ ScaleSpace<I> *ss = new ScaleSpace<I>(octave_level, blur_level);
+ DoGSpace<I, S>* dogs = new DoGSpace<I, S>(ss);
+ std::list<Keypoint> keypoints;
+
I original;
- C extrema, improved;
- bool black = true;
+ C extrema, extrema_b,
+ extrema2, extrema2_b;
+ //improved, improved_b;
+
+ io::pgm::load(original, source.c_str());
- io::pgm::load(original, "images/keith.pbm");
+ initialize(extrema_b, original);
+ initialize(extrema2_b, original);
+ //initialize(improved_b, original);
- if (black)
- {
- initialize(extrema, original);
- initialize(improved, original);
- }
- else
- {
- extrema = data::convert(value::rgb8(), original);
- improved = data::convert(value::rgb8(), original);
- }
+ extrema = data::convert(value::rgb8(), original);
+ extrema2 = data::convert(value::rgb8(), original);
+ //improved = data::convert(value::rgb8(), original);
// Localization
- buildScaleSpace(scaleSpace, original, octave_level, blur_level);
- buildDifferenceOfGaussianSpace(scaleSpace, dogSpace);
- buildExtrema(extrema, dogSpace, keypoints);
- discardLowContrastKeypoints(dogSpace, keypoints, improved);
+ ss->build(original);
+ //ss->save();
+ dogs->build();
+ //dogs->save();
+ dogs->findKeypoints(keypoints);
+
+ writeKeypoints(keypoints, extrema, extrema_b);
+
+ dogs->discardLowContrastKeypoints(keypoints);
+ dogs->eliminateEdgeResponses(keypoints);
- // Processing
- // TODO
+ writeKeypoints(keypoints, extrema2, extrema2_b);
+ // Save
io::ppm::save(extrema, "output/extrema.ppm");
- io::ppm::save(improved, "output/extrema_improved.ppm");
+ io::ppm::save(extrema_b, "output/extrema_b.ppm");
+ io::ppm::save(extrema2, "output/extrema2.ppm");
+ io::ppm::save(extrema2_b, "output/extrema2_b.ppm");
+ //io::ppm::save(improved, "output/extrema_improved.ppm");
+ //io::ppm::save(improved_b, "output/extrema_b_improved.ppm");
return 0;
}
diff --git a/milena/sandbox/anthony/src/matrix.hh b/milena/sandbox/anthony/src/matrix.hh
index de83912..c42de59 100644
--- a/milena/sandbox/anthony/src/matrix.hh
+++ b/milena/sandbox/anthony/src/matrix.hh
@@ -16,6 +16,7 @@ class Matrix
bool inverse(Matrix& inversed);
void absolute();
bool isMajorOffset();
+ float determinant();
// Operators
Matrix operator*(Matrix& right);
@@ -31,8 +32,6 @@ class Matrix
std::vector< std::vector<float> > m;
unsigned height;
unsigned width;
-
- float determinant();
};
#endif /* ! MATRIX_HH */
diff --git a/milena/sandbox/anthony/new-implem/src/scale-space.hh b/milena/sandbox/anthony/src/scale-space.hh
similarity index 100%
rename from milena/sandbox/anthony/new-implem/src/scale-space.hh
rename to milena/sandbox/anthony/src/scale-space.hh
diff --git a/milena/sandbox/anthony/new-implem/src/vrac.hh b/milena/sandbox/anthony/src/vrac.hh
similarity index 100%
rename from milena/sandbox/anthony/new-implem/src/vrac.hh
rename to milena/sandbox/anthony/src/vrac.hh
--
1.7.10.4
