[Olena-patches] olena: olena-2.0-595-gc54e2d1 Improve extremum

--- milena/sandbox/anthony/Makefile | 2 +- milena/sandbox/anthony/src/keypoint.cc | 39 +++++- milena/sandbox/anthony/src/keypoint.hh | 29 +++-- milena/sandbox/anthony/src/main.cc | 226 ++++++++++++++++++++++++------- milena/sandbox/anthony/src/matrix.cc | 30 ++++- milena/sandbox/anthony/src/matrix.hh | 5 +- 6 files changed, 259 insertions(+), 72 deletions(-) diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/Makefile index 83ee52b..2953013 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/matrix.cc src/keypoint.cc src/main.cc +SRC=src/main.cc src/matrix.cc src/keypoint.cc OUTPUT=a.out all: scale diff --git a/milena/sandbox/anthony/src/keypoint.cc b/milena/sandbox/anthony/src/keypoint.cc index bf0a5f5..937e78d 100644 --- a/milena/sandbox/anthony/src/keypoint.cc +++ b/milena/sandbox/anthony/src/keypoint.cc @@ -1,8 +1,8 @@ #include "keypoint.hh" // Constructor -Keypoint::Keypoint(unsigned i, unsigned j, - unsigned scale, unsigned size, bool maximum) +Keypoint::Keypoint(int i, int j, + int scale, int size, bool maximum) { this->i = i; this->j = j; @@ -11,8 +11,35 @@ Keypoint::Keypoint(unsigned i, unsigned j, 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 -unsigned Keypoint::getI() { return i; } -unsigned Keypoint::getJ() { return j; } -unsigned Keypoint::getScale() { return scale; } -unsigned Keypoint::getSize() { return size; } +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 798a980..af24a4c 100644 --- a/milena/sandbox/anthony/src/keypoint.hh +++ b/milena/sandbox/anthony/src/keypoint.hh @@ -1,26 +1,35 @@ #ifndef KEYPOINT_HH # define KEYPOINT_HH +//# include <mln/core/alias/point2d.hh> # include <ostream> +# include "matrix.hh" class Keypoint { public: // Constructor - Keypoint(unsigned i, unsigned j, - unsigned scale, unsigned size, bool maximum); + Keypoint(int i, int j, + int scale, int size, bool maximum); + + bool add(Matrix& k, Keypoint& old); + //mln::point2d getPoint(); // Getters - unsigned getI(); - unsigned getJ(); - unsigned getScale(); - unsigned getSize(); + int getI(); + int getJ(); + int getScale(); + int getSize(); + + void setI(int offsetI); + void setJ(int offsetJ); + void setScale(int offsetScale); private: - unsigned i; - unsigned j; - unsigned scale; - unsigned size; + int i; + int j; + int scale; + int size; bool maximum; }; diff --git a/milena/sandbox/anthony/src/main.cc b/milena/sandbox/anthony/src/main.cc index 9cc956b..fae3c64 100644 --- a/milena/sandbox/anthony/src/main.cc +++ b/milena/sandbox/anthony/src/main.cc @@ -472,18 +472,11 @@ bool hasEdgeResponse(const I& dog, const P& p, const T& center) // Build the extrema image template<typename I, typename C> -void buildExtrema(C extrema, - I original, +void buildExtrema(C& extrema, + const I& original, std::vector< std::vector<I> >& dogSpace, std::vector<Keypoint>& keypoints) { - const bool black = false; - - if (black) - initialize(extrema, original); - else - extrema = data::convert(value::rgb8(), original); - // Number of scales (gaussian meaning) unsigned scales = dogSpace.at(0).size(); @@ -523,7 +516,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; + extrema(pOriginal) = literal::red; keypoints.push_back(Keypoint(p.row(), p.col(), i, j, false)); } } @@ -543,75 +536,191 @@ void buildExtrema(C extrema, } } } - - io::ppm::save(extrema, "output/extrema.ppm"); } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - +// Compute the first order matrix for interpolation template<typename I> -Matrix *computeFirstOrderMatrix(const std::vector< std::vector<I> >& dogSpace, - const Keypoint& k) +void computeFirstOrderMatrix(const std::vector< std::vector<I> >& dogSpace, + Keypoint& k, + Matrix& matrix) { - Matrix *matrix = new Matrix(3, 1); + 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); - return matrix; + // 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> -Matrix *computeSecondOrderMatrix(const std::vector< std::vector<I> >& dogSpace, - const Keypoint& k) +void computeSecondOrderMatrix(const std::vector< std::vector<I> >& dogSpace, + Keypoint& k, + Matrix& matrix) { - Matrix *matrix = new Matrix(3, 3); + matrix.initialize(); - return matrix; -} + 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); -bool isBetter(Matrix offset) -{ - return (offset.get(0, 0) > 0.5 - || offset.get(0, 1) > 0.5 - || offset.get(0, 2) > 0.5); -} + /* Point naming convention + * n5 n1 n6 + * n2 c n3 + * n7 n4 n8 + */ + point2d c(k.getI(), k.getJ()); -void changeKeypoint(Keypoint& k, Matrix& offset) -{ + 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) -// FIXME Move to object representation (size-templated matrix maybe ?) -template<typename I> +template<typename I, typename C> void discardLowContrastKeypoints(const std::vector< std::vector<I> >& dogSpace, - std::vector<Keypoint>& keypoints) + std::vector<Keypoint>& keypoints, + C& extrema) { + Keypoint k_old(0, 0, 0, 0, true); + 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); + Matrix firstOrder(3, 1); + Matrix secondOrder(3, 3); + Matrix inversed(secondOrder.getHeight(), secondOrder.getWidth()); - if (inversible) + computeFirstOrderMatrix(dogSpace, k, firstOrder); + computeSecondOrderMatrix(dogSpace, k, secondOrder); + bool hasInversed = secondOrder.inverse(inversed); + + if (hasInversed) { - Matrix offset = (*fom) * (*som); - - if (isBetter(offset)) + Matrix x = inversed * firstOrder; + bool isMajorOffset = x.isMajorOffset(); + + if (isMajorOffset) + { + std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + std::cout << "First order matrix:" << std::endl; + firstOrder.print(); + + std::cout << "Second order matrix:" << std::endl; + secondOrder.print(); + + std::cout << "Offset: " << std::endl; + x.print(); + + std::cout << "k_old: " << k_old.getI() << " " << k_old.getJ() << std::endl; + std::cout << "k: " << k.getI() << " " << k.getJ() << std::endl; + bool hasMoved = k.add(x, k_old); + std::cout << "k++: " << k.getI() << " " << k.getJ() + << "\thasMoved: " << hasMoved << std::endl; + + point2d p(k.getI(), k.getJ()); + std::cout << "Final point: " << extrema.domain() << "\t" << p << std::endl; + + if (hasMoved && extrema.has(p) + && p.row() >= 0 && p.row() < geom::max_row(extrema) + && p.col() >= 0 && p.col() < geom::max_col(extrema)) + { + std::cout << "VALID" << std::endl; + extrema(p) = literal::red; + k_old = keypoints.at(i); + keypoints.at(i--) = k; + } + + std::cout << std::endl; + } + else { - changeKeypoint(k, offset); + point2d p(k.getI(), k.getJ()); + extrema(p) = literal::red; } } } } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - // Main entry point int main(int argc, char** argv) { @@ -624,18 +733,33 @@ int main(int argc, char** argv) std::vector< std::vector<I> > dogSpace; std::vector<Keypoint> keypoints; I original; - C extrema; + C extrema, improved; + bool black = true; io::pgm::load(original, "images/lena.pgm"); + if (black) + { + initialize(extrema, original); + initialize(improved, original); + } + else + { + extrema = data::convert(value::rgb8(), original); + improved = data::convert(value::rgb8(), original); + } + // Localization buildScaleSpace(scaleSpace, original, octave_level, blur_level); buildDifferenceOfGaussianSpace(scaleSpace, dogSpace); buildExtrema(extrema, original, dogSpace, keypoints); - //discardLowContrastKeypoints(dogSpace, keypoints); + discardLowContrastKeypoints(dogSpace, keypoints, improved); // Processing // TODO + io::ppm::save(extrema, "output/extrema.ppm"); + io::ppm::save(improved, "output/extrema_improved.ppm"); + return 0; } diff --git a/milena/sandbox/anthony/src/matrix.cc b/milena/sandbox/anthony/src/matrix.cc index 1c5784d..ee00ffc 100644 --- a/milena/sandbox/anthony/src/matrix.cc +++ b/milena/sandbox/anthony/src/matrix.cc @@ -46,7 +46,6 @@ void Matrix::print() std::cout << std::endl; } - float Matrix::determinant() { float d1 = m.at(0).at(0) * m.at(1).at(1) * m.at(2).at(2); @@ -99,7 +98,6 @@ bool Matrix::inverse(Matrix& inversed) for (unsigned j = 0; j < comatrix.getHeight(); ++j) comatrix.set(i, j, comatrix.get(i,j) / det); - comatrix.print(); inversed = comatrix; canInverse = true; @@ -120,7 +118,33 @@ Matrix Matrix::operator*(Matrix& right) 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; } diff --git a/milena/sandbox/anthony/src/matrix.hh b/milena/sandbox/anthony/src/matrix.hh index b2ca124..de83912 100644 --- a/milena/sandbox/anthony/src/matrix.hh +++ b/milena/sandbox/anthony/src/matrix.hh @@ -13,8 +13,9 @@ class Matrix void initialize(); void initialize(float *values); void print(); - float determinant(); bool inverse(Matrix& inversed); + void absolute(); + bool isMajorOffset(); // Operators Matrix operator*(Matrix& right); @@ -30,6 +31,8 @@ class Matrix std::vector< std::vector<float> > m; unsigned height; unsigned width; + + float determinant(); }; #endif /* ! MATRIX_HH */ -- 1.7.2.5