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

28 May 2013

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

    

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

Anthony Seure