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