[Olena-patches] olena: olena-2.0-599-g533f3a6 [SIFT] Replace old implementation

--- 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