[Olena-patches] olena: olena-2.0-594-g3748427 Add matrix class to handle Taylor approximation (dirty version)

--- milena/sandbox/anthony/Makefile | 2 +- milena/sandbox/anthony/src/main.cc | 78 ++++++++++++++++++--- milena/sandbox/anthony/src/matrix.cc | 130 ++++++++++++++++++++++++++++++++++ milena/sandbox/anthony/src/matrix.hh | 35 +++++++++ 4 files changed, 234 insertions(+), 11 deletions(-) create mode 100644 milena/sandbox/anthony/src/matrix.cc create mode 100644 milena/sandbox/anthony/src/matrix.hh diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/Makefile index ff62f9e..83ee52b 100644 --- a/milena/sandbox/anthony/Makefile +++ b/milena/sandbox/anthony/Makefile @@ -4,7 +4,7 @@ CFLAGS=-Wall -Werror CLIBS=-I../../ CLEAN=*.o output/* log -SRC=src/main.cc src/keypoint.cc +SRC=src/matrix.cc src/keypoint.cc src/main.cc OUTPUT=a.out all: scale diff --git a/milena/sandbox/anthony/src/main.cc b/milena/sandbox/anthony/src/main.cc index eea1c86..9cc956b 100644 --- a/milena/sandbox/anthony/src/main.cc +++ b/milena/sandbox/anthony/src/main.cc @@ -9,6 +9,7 @@ #include <cmath> #include "keypoint.hh" +#include "matrix.hh" using namespace mln; @@ -523,7 +524,7 @@ void buildExtrema(C extrema, mln_psite(I) pOriginal(p.row() * pow(2, (j-1)), p.col() * pow(2, (j-1))); extrema(pOriginal) = literal::green; - keypoints.push_back(Keypoint(p.row(), p.col(), j, i, false)); + keypoints.push_back(Keypoint(p.row(), p.col(), i, j, false)); } } else if (center >= max) @@ -536,7 +537,7 @@ void buildExtrema(C extrema, 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(), j, i, true)); + keypoints.push_back(Keypoint(p.row(), p.col(), i, j, true)); } } } @@ -546,13 +547,71 @@ void buildExtrema(C extrema, io::ppm::save(extrema, "output/extrema.ppm"); } +/*****************************************************************************/ +/*****************************************************************************/ +/*****************************************************************************/ + template<typename I> -void discardEdgeResponses(std::vector< std::vector<I> >& dogSpace, - std::vector<Keypoint>& keypoints) +Matrix *computeFirstOrderMatrix(const std::vector< std::vector<I> >& dogSpace, + const Keypoint& k) +{ + Matrix *matrix = new Matrix(3, 1); + + return matrix; +} + +template<typename I> +Matrix *computeSecondOrderMatrix(const std::vector< std::vector<I> >& dogSpace, + const Keypoint& k) +{ + Matrix *matrix = new Matrix(3, 3); + + return matrix; +} + +bool isBetter(Matrix offset) { - // FIXME Move work from buildExtrema here + return (offset.get(0, 0) > 0.5 + || offset.get(0, 1) > 0.5 + || offset.get(0, 2) > 0.5); } +void changeKeypoint(Keypoint& k, Matrix& offset) +{ +} + +// Discard low contrast keypoints thanks to taylor interpolation +// See Brown and Lowe paper (2002) +// FIXME Move to object representation (size-templated matrix maybe ?) +template<typename I> +void discardLowContrastKeypoints(const std::vector< std::vector<I> >& dogSpace, + std::vector<Keypoint>& keypoints) +{ + for (unsigned i = 0; i < keypoints.size(); ++i) + { + Keypoint k = keypoints.at(i); + + Matrix *fom = computeFirstOrderMatrix(dogSpace, k); + Matrix *som = computeSecondOrderMatrix(dogSpace, k); + Matrix inversed(som->getHeight(), som->getWidth()); + bool inversible = som->inverse(inversed); + + if (inversible) + { + Matrix offset = (*fom) * (*som); + + if (isBetter(offset)) + { + changeKeypoint(k, offset); + } + } + } +} + +/*****************************************************************************/ +/*****************************************************************************/ +/*****************************************************************************/ + // Main entry point int main(int argc, char** argv) { @@ -569,15 +628,14 @@ int main(int argc, char** argv) io::pgm::load(original, "images/lena.pgm"); + // Localization buildScaleSpace(scaleSpace, original, octave_level, blur_level); buildDifferenceOfGaussianSpace(scaleSpace, dogSpace); buildExtrema(extrema, original, dogSpace, keypoints); - // TODO Interpolation of nearby data for accurate position - // TODO Discard low-contrast keypoints - //discardEdgeResponses(dogSpace, keypoints); - - std::cout << keypoints.size() << std::endl; + //discardLowContrastKeypoints(dogSpace, keypoints); + // Processing + // TODO return 0; } diff --git a/milena/sandbox/anthony/src/matrix.cc b/milena/sandbox/anthony/src/matrix.cc new file mode 100644 index 0000000..1c5784d --- /dev/null +++ b/milena/sandbox/anthony/src/matrix.cc @@ -0,0 +1,130 @@ +#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); + + comatrix.print(); + 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; +}; + +// 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/src/matrix.hh b/milena/sandbox/anthony/src/matrix.hh new file mode 100644 index 0000000..b2ca124 --- /dev/null +++ b/milena/sandbox/anthony/src/matrix.hh @@ -0,0 +1,35 @@ +#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(); + float determinant(); + bool inverse(Matrix& inversed); + + // 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 */ -- 1.7.2.5