---
milena/sandbox/anthony/Makefile | 2 +-
milena/sandbox/anthony/opencv/Makefile | 5 +-
milena/sandbox/anthony/opencv/main.cpp | 32 +-
milena/sandbox/anthony/src/dog-space.hh | 690 +++++++++++++++--------------
milena/sandbox/anthony/src/main.cc | 48 +-
milena/sandbox/anthony/src/scale-space.hh | 21 +-
6 files changed, 412 insertions(+), 386 deletions(-)
diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/Makefile
index 44edf3b..edd3807 100644
--- a/milena/sandbox/anthony/Makefile
+++ b/milena/sandbox/anthony/Makefile
@@ -1,4 +1,4 @@
-CCACHE=
+CCACHE=ccache
CC=g++
CFLAGS=-Wall -Wextra
CLIBS=-I../../
diff --git a/milena/sandbox/anthony/opencv/Makefile b/milena/sandbox/anthony/opencv/Makefile
index 8e438f0..d9a8be4 100644
--- a/milena/sandbox/anthony/opencv/Makefile
+++ b/milena/sandbox/anthony/opencv/Makefile
@@ -2,7 +2,7 @@ CCACHE=
CXX=clang++
SRC=main.cpp
LIBS=-L/usr/local/include -lopencv_highgui -lopencv_core -lopencv_features2d -lopencv_nonfree -D_GNU_SOURCE
-OUTPUT=cv
+OUTPUT=a.out
CLEAN=${OUTPUT} *.jpg
all: cv
@@ -10,6 +10,9 @@ all: cv
cv:
${CCACHE} ${CXX} -O3 ${LIBS} ${SRC} -o ${OUTPUT}
+pyramid:
+ ${CCACHE} ${CXX} -O3 ${LIBS} pyramid.cpp -o ${OUTPUT}
+
debug:
${CCACHE} ${CXX} -g -ggdb -O0 ${LIBS} ${SRC} -o ${OUTPUT}
diff --git a/milena/sandbox/anthony/opencv/main.cpp b/milena/sandbox/anthony/opencv/main.cpp
index 2de37c1..a7e5713 100644
--- a/milena/sandbox/anthony/opencv/main.cpp
+++ b/milena/sandbox/anthony/opencv/main.cpp
@@ -10,14 +10,30 @@
// MAIN ENTRY POINT
int main(int argc, char* argv[])
{
- const cv::Mat input = cv::imread("../images/keith.pbm");
- cv::SiftFeatureDetector detector;
- std::vector<cv::KeyPoint> keypoints;
- detector.detect(input, keypoints);
-
- cv::Mat output;
- cv::drawKeypoints(input, keypoints, output);
- cv::imwrite("sift.jpg", output);
+ const cv::Mat input1 = cv::imread("../images/keith.pbm");
+ const cv::Mat input2 = cv::imread("../images/tshirt.pgm");
+
+ {
+ cv::SiftFeatureDetector detector;
+ std::vector<cv::KeyPoint> keypoints;
+ detector.detect(input1, keypoints);
+ std::cout << keypoints.size() << std::endl;
+
+ cv::Mat output;
+ cv::drawKeypoints(input1, keypoints, output);
+ cv::imwrite("sift1.jpg", output);
+ }
+
+ {
+ cv::SiftFeatureDetector detector;
+ std::vector<cv::KeyPoint> keypoints;
+ detector.detect(input2, keypoints);
+ std::cout << keypoints.size() << std::endl;
+
+ cv::Mat output;
+ cv::drawKeypoints(input2, keypoints, output);
+ cv::imwrite("sift2.jpg", output);
+ }
return 0;
diff --git a/milena/sandbox/anthony/src/dog-space.hh b/milena/sandbox/anthony/src/dog-space.hh
index da35820..9c4fee9 100644
--- a/milena/sandbox/anthony/src/dog-space.hh
+++ b/milena/sandbox/anthony/src/dog-space.hh
@@ -12,6 +12,7 @@
#include <mln/core/alias/neighb2d.hh>
#include <mln/literal/all.hh>
#include <mln/value/int_s16.hh>
+#include <mln/value/int_u8.hh>
#include "scale-space.hh"
#include "keypoint.hh"
@@ -22,377 +23,384 @@ using namespace mln;
template<typename I, typename S>
class DoGSpace : protected ScaleSpace<S>
{
- public:
- DoGSpace(ScaleSpace<I>* _ss)
+ public:
+ DoGSpace(ScaleSpace<I>* ss)
+ : ss_(ss)
+ { /* NOTHING */ }
+
+ ~DoGSpace()
+ {
+ ScaleSpace<S>::~ScaleSpace();
+ }
+
+ void build(void)
+ {
+ ScaleSpace<S>::o = ss_->octave_count();
+ ScaleSpace<S>::g = ss_->gradient_count() - 1;
+
+ typedef typename std::vector< std::vector<I> >::iterator PIT;
+ typedef typename std::vector<I>::iterator IIT;
+
+ std::vector< std::vector<I> > sspyramid = ss_->getPyramid();
+ S dog;
+
+ for (PIT pit = sspyramid.begin(); pit != sspyramid.end(); ++pit)
+ {
+ std::vector<S> octave;
+ for (IIT iit = pit->begin(); (iit + 1) != pit->end(); ++iit)
{
- ss = _ss;
- }
+ initialize(dog, *iit);
+ mln_piter(I) p(iit->domain());
- ~DoGSpace()
- {
- ScaleSpace<S>::~ScaleSpace();
- }
+ for_all(p)
+ {
+ value::int_u8 up = (*iit)(p);
+ value::int_u8 down = (*(iit+1))(p);
- void build(void)
+ dog(p) = static_cast<value::int_s16>(up) - static_cast<value::int_s16>(down);
+ }
+ octave.push_back(dog);
+ }
+ ScaleSpace<S>::pyramid.push_back(octave);
+ }
+ }
+
+ void findKeypoints(std::list<Keypoint>& keypoints)
+ {
+ typedef typename std::vector< std::vector<S> >::iterator PIT;
+ typedef typename std::vector<S>::iterator IIT;
+
+ unsigned octave = 0;
+ for (PIT pit = ScaleSpace<S>::pyramid.begin(); pit != ScaleSpace<S>::pyramid.end(); ++pit)
+ {
+ unsigned gradient = 1;
+ for (IIT iit = pit->begin() + 1; (iit + 1) != pit->end(); ++iit)
{
- ScaleSpace<S>::o = ss->octave_count();
- ScaleSpace<S>::g = ss->gradient_count() - 1;
+ S upper = *(iit-1);
+ S current = *iit;
+ S lower = *(iit+1);
- typedef typename std::vector< std::vector<I> >::iterator PIT;
- typedef typename std::vector<I>::iterator IIT;
+ mln_piter(I) p(current.domain());
- std::vector< std::vector<I> > sspyramid = ss->getPyramid();
- S dog;
+ for_all(p)
+ {
+ value::int_s16 center = current(p);
+ value::int_s16 min = 255;
+ value::int_s16 max = -256;
+ mln_niter(neighb2d) n(c8(), p);
- for (PIT pit = sspyramid.begin(); pit != sspyramid.end(); ++pit)
+ int extrema = isExtrema(current, center, n);
+
+ if (extrema == -1) // Is a minimum
{
- std::vector<S> octave;
- for (IIT iit = pit->begin(); (iit + 1) != pit->end(); ++iit)
- {
- initialize(dog, *iit);
- mln_piter(I) p(iit->domain());
-
- for_all(p)
- {
- value::int_s16 up = (value::int_s16) (*iit)(p);
- value::int_s16 down = (value::int_s16) (*(iit+1))(p);
- value::int_s16 diff = down - up;
- //value::int_s16 diff = ( (int) (*iit)(p) ) - ( (int) ((*(iit+1))(p)) );
-
- //std::cout << up << " - " << down << " = " << diff << std::endl;
-
- dog(p) = diff;
- }
- octave.push_back(dog);
- }
- ScaleSpace<S>::pyramid.push_back(octave);
+ findMin(min, lower, upper, p, n);
+ if (center < min)
+ keypoints.push_back(Keypoint(p.col(), p.row(), octave, gradient, false));
}
- }
-
- void findKeypoints(std::list<Keypoint>& keypoints)
- {
- typedef typename std::vector< std::vector<S> >::iterator PIT;
- typedef typename std::vector<S>::iterator IIT;
-
- unsigned octave = 0;
- for (PIT pit = ScaleSpace<S>::pyramid.begin(); pit != ScaleSpace<S>::pyramid.end(); ++pit)
+ else if (extrema == 1) // Is a maximum
{
- unsigned gradient = 1;
- for (IIT iit = pit->begin() + 1; (iit + 1) != pit->end(); ++iit)
- {
- S upper = *(iit-1);
- S current = *iit;
- S lower = *(iit+1);
-
- mln_piter(I) p(current.domain());
-
- for_all(p)
- {
- value::int_s16 center = current(p);
- value::int_s16 min = 255;
- value::int_s16 max = -256;
- mln_niter(neighb2d) n(c8(), p);
-
- int extrema = isExtrema(current, center, n);
-
- if (extrema == -1) // Is a minimum
- {
- findMin(min, lower, upper, p, n);
- if (center < min)
- keypoints.push_back(Keypoint(p.col(), p.row(), octave, gradient, false));
- }
- else if (extrema == 1) // Is a maximum
- {
- findMax(max, lower, upper, p, n);
- if (center > max)
- keypoints.push_back(Keypoint(p.col(), p.row(), octave, gradient, true));
- }
- }
-
- ++gradient;
- }
-
- ++octave;
+ findMax(max, lower, upper, p, n);
+ if (center > max)
+ keypoints.push_back(Keypoint(p.col(), p.row(), octave, gradient, true));
}
+ }
+
+ ++gradient;
}
- void discardLowContrastKeypoints(std::list<Keypoint>& keypoints)
+ ++octave;
+ }
+ }
+
+ void discardLowContrastKeypoints(std::list<Keypoint>& keypoints)
+ {
+ const static unsigned max_interpolation_depth = 5;
+ unsigned count = 0;
+ bool enough = false;
+ std::list<Keypoint>::iterator k = keypoints.begin();
+ std::list<Keypoint>::iterator end = keypoints.end();
+ Matrix* offset = NULL;
+ Matrix* previous_offset = NULL;
+
+ while (k != end)
+ {
+ while (!enough && count < max_interpolation_depth)
{
- const static unsigned max_interpolation_depth = 5;
- unsigned count = 0;
- bool enough = false;
- std::list<Keypoint>::iterator k = keypoints.begin();
- Matrix* offset = NULL;
- Matrix* previous_offset = NULL;
-
- while (k != keypoints.end())
+ unsigned octave = k->getOctave();
+ unsigned gradient = k->getGradient();
+ S dog = ScaleSpace<S>::pyramid[octave][gradient];
+ S dog_up = ScaleSpace<S>::pyramid[octave][gradient + 1];
+ S dog_down = ScaleSpace<S>::pyramid[octave][gradient - 1];
+ int col = k->getX();
+ int row = k->getY();
+
+ if (row < 0 || col < 0
+ || static_cast<unsigned>(row) >= dog.nrows()
+ || static_cast<unsigned>(col) >= dog.ncols())
+ {
+ count = max_interpolation_depth;
+ break;
+ }
+
+ /* Set the neighbors of the current level:
+ * n5 n1 n6
+ * n2 p n3
+ * n7 n4 n8
+ */
+ float p = dog(mln_psite(S)(row, col));
+
+ float n1 = dog(mln_psite(S)(row - 1, col));
+ float n2 = dog(mln_psite(S)(row, col - 1));
+ float n3 = dog(mln_psite(S)(row, col + 1));
+ float n4 = dog(mln_psite(S)(row + 1, col));
+
+ float n5 = dog(mln_psite(S)(row - 1, col - 1));
+ float n6 = dog(mln_psite(S)(row - 1, col + 1));
+ float n7 = dog(mln_psite(S)(row + 1, col - 1));
+ float n8 = dog(mln_psite(S)(row + 1, col + 1));
+
+ /* Set the neighbors of the upper level:
+ * u1
+ * u2 u u3
+ * u4
+ */
+ float u = dog_up(mln_psite(S)(row, col));
+ float u1 = dog_up(mln_psite(S)(row - 1, col));
+ float u2 = dog_up(mln_psite(S)(row, col - 1));
+ float u3 = dog_up(mln_psite(S)(row, col + 1));
+ float u4 = dog_up(mln_psite(S)(row + 1, col));
+
+ /* Set the neighbors of the lower level:
+ * l1
+ * l2 l l3
+ * l4
+ */
+ float l = dog_down(mln_psite(S)(row, col));
+ float l1 = dog_down(mln_psite(S)(row - 1, col));
+ float l2 = dog_down(mln_psite(S)(row, col - 1));
+ float l3 = dog_down(mln_psite(S)(row, col + 1));
+ float l4 = dog_down(mln_psite(S)(row + 1, col));
+
+ // Compute the offset
+
+ float dxx = n3 - 2.0 * p + n2;
+ float dyy = n4 - 2.0 * p + n1;
+ float dss = u - 2.0 * p + l;
+
+ float dxy = (n8 - n6 - n7 + n5) / 4.0;
+ float dxs = (u3 - u2 - l3 + l2) / 4.0;
+ float dys = (u4 - u1 - l4 + l1) / 4.0;
+
+ Matrix m1(3, 1);
+ m1.initialize();
+ m1.set(0, 0, n3 - n2);
+ m1.set(1, 0, n4 - n1);
+ m1.set(2, 0, u - l);
+
+ Matrix m2(3, 3);
+ m2.initialize();
+ m2.set(0, 0, dxx); m2.set(0, 1, dxy); m2.set(0, 2, dxs);
+ m2.set(1, 0, dxy); m2.set(1, 1, dyy); m2.set(1, 2, dys);
+ m2.set(2, 0, dxs); m2.set(2, 1, dys); m2.set(2, 2, dss);
+
+ previous_offset = offset;
+ offset = m2.solve3x3(m1);
+
+ if (offset == NULL)
+ {
+ offset = previous_offset;
+ enough = true;
+ }
+ else
+ {
+ if (std::abs(offset->get(0,0)) > 0.5 ||
+ std::abs(offset->get(1,0)) > 0.5 ||
+ std::abs(offset->get(2,0)) > 0.5)
{
- while (!enough && count < max_interpolation_depth)
- {
- std::cout << "test1" << std::endl;
- unsigned octave = k->getOctave();
- unsigned gradient = k->getGradient();
- S dog = ScaleSpace<S>::pyramid[octave][gradient];
- S dog_up = ScaleSpace<S>::pyramid[octave][gradient + 1];
- S dog_down = ScaleSpace<S>::pyramid[octave][gradient - 1];
- unsigned col = k->getX();
- unsigned row = k->getY();
- std::cout << "test2" << std::endl;
-
- /* Set the neighbors of the current level:
- * n5 n1 n6
- * n2 p n3
- * n7 n4 n8
- */
- float p = dog(mln_psite(S)(row, col));
-
- float n1 = dog(mln_psite(S)(row - 1, col));
- float n2 = dog(mln_psite(S)(row, col - 1));
- float n3 = dog(mln_psite(S)(row, col + 1));
- float n4 = dog(mln_psite(S)(row + 1, col));
-
- float n5 = dog(mln_psite(S)(row - 1, col - 1));
- float n6 = dog(mln_psite(S)(row - 1, col + 1));
- float n7 = dog(mln_psite(S)(row + 1, col - 1));
- float n8 = dog(mln_psite(S)(row + 1, col + 1));
-
- /* Set the neighbors of the upper level:
- * u1
- * u2 u u3
- * u4
- */
- float u = dog_up(mln_psite(S)(row, col));
- float u1 = dog_up(mln_psite(S)(row - 1, col));
- float u2 = dog_up(mln_psite(S)(row, col - 1));
- float u3 = dog_up(mln_psite(S)(row, col + 1));
- float u4 = dog_up(mln_psite(S)(row + 1, col));
-
- /* Set the neighbors of the lower level:
- * l1
- * l2 l l3
- * l4
- */
- float l = dog_down(mln_psite(S)(row, col));
- float l1 = dog_down(mln_psite(S)(row - 1, col));
- float l2 = dog_down(mln_psite(S)(row, col - 1));
- float l3 = dog_down(mln_psite(S)(row, col + 1));
- float l4 = dog_down(mln_psite(S)(row + 1, col));
-
- // Compute the offset
-
- float dxx = n3 - 2.0 * p + n2;
- float dyy = n4 - 2.0 * p + n1;
- float dss = u - 2.0 * p + l;
-
- float dxy = (n8 - n6 - n7 + n5) / 4.0;
- float dxs = (u3 - u2 - l3 + l2) / 4.0;
- float dys = (u4 - u1 - l4 + l1) / 4.0;
-
- std::cout << "test3" << std::endl;
- Matrix m1(3, 1);
- m1.initialize();
- m1.set(0, 0, n3 - n2);
- m1.set(1, 0, n4 - n1);
- m1.set(2, 0, u - l);
-
- Matrix m2(3, 3);
- m2.initialize();
- m2.set(0, 0, dxx); m2.set(0, 1, dxy); m2.set(0, 2, dxs);
- m2.set(1, 0, dxy); m2.set(1, 1, dyy); m2.set(1, 2, dys);
- m2.set(2, 0, dxs); m2.set(2, 1, dys); m2.set(2, 2, dss);
-
- std::cout << "test4" << std::endl;
- previous_offset = offset;
- offset = m2.solve3x3(m1);
- std::cout << "test5" << std::endl;
-
- if (offset == NULL)
- {
- offset = previous_offset;
- enough = true;
- }
- else
- {
- if (std::abs(offset->get(0,0)) > 0.5 ||
- std::abs(offset->get(1,0)) > 0.5 ||
- std::abs(offset->get(2,0)) > 0.5)
- {
- std::cout << "test51" << std::endl;
- ++count;
- k->setX(k->getX() + (int) offset->get(0, 0));
- k->setY(k->getX() + (int) offset->get(1, 0));
- k->setGradient(k->getGradient() + (int) offset->get(1, 0));
-
- //FIXME Check boundaries
-
- if (k->getGradient() < 1 ||
- k->getGradient() >= ScaleSpace<S>::pyramid[0].size() - 1)
- count = max_interpolation_depth; // Remove keypoint
- }
- else
- {
- std::cout << "test52" << std::endl;
- enough = true;
- }
- }
- }
-
- std::cout << "test6" << std::endl;
- if (count >= max_interpolation_depth)
- {
- std::cout << "test61" << std::endl;
- k = keypoints.erase(k);
- }
- else
- {
- std::cout << "test62" << std::endl;
- unsigned octave = k->getOctave();
- unsigned gradient = k->getGradient();
-
- S dog = ScaleSpace<S>::pyramid[octave][gradient];
- S dog_up = ScaleSpace<S>::pyramid[octave][gradient + 1];
- S dog_down = ScaleSpace<S>::pyramid[octave][gradient - 1];
- unsigned col = k->getX();
- unsigned row = k->getY();
-
- /* Set the neighbors of the current level:
- * n1
- * n2 p n3
- * n4
- */
- float p = dog(mln_psite(S)(row, col));
-
- float n1 = dog(mln_psite(S)(row - 1, col));
- float n2 = dog(mln_psite(S)(row, col - 1));
- float n3 = dog(mln_psite(S)(row, col + 1));
- float n4 = dog(mln_psite(S)(row + 1, col));
-
- // Set the neighbors of the upper/lower level:
- float u = dog_up(mln_psite(S)(row, col));
- float l = dog_down(mln_psite(S)(row, col));
-
- Matrix firstOrder(3, 1);
- firstOrder.initialize();
-
- /************ WARNING ************/
- /* Maybe missing /255 */
- /*********************************/
-
- firstOrder.set(0, 0, n3 - n2);
- firstOrder.set(1, 0, n4 - n1);
- firstOrder.set(2, 0, u - l);
-
- float product = firstOrder.dotProduct(offset);
- float contrast = p + product / 2.0;
-
- if (contrast < 0.03)
- k = keypoints.erase(k);
- else
- k++;
- }
-
- enough = false;
- count = 0;
- }
- }
+ ++count;
+ k->setX(k->getX() + (int) offset->get(0, 0));
+ k->setY(k->getX() + (int) offset->get(1, 0));
+ k->setGradient(k->getGradient() + (int) offset->get(1, 0));
- void eliminateEdgeResponses(std::list<Keypoint>& keypoints)
- {
- const static float bound_ratio = 12.1;
+ //FIXME Check boundaries
- for (std::list<Keypoint>::iterator k = keypoints.begin();
- k != keypoints.end(); ++k)
+ if (k->getGradient() < 1 ||
+ k->getGradient() >= ScaleSpace<S>::pyramid[0].size() - 1)
+ count = max_interpolation_depth; // Remove keypoint
+ }
+ else
{
- unsigned octave = k->getOctave();
- unsigned gradient = k->getGradient();
- S dog = ScaleSpace<S>::pyramid[octave][gradient];
- unsigned col = k->getX();
- unsigned row = k->getY();
-
- /* Set the neighbors:
- * n5 n1 n6
- * n2 p n3
- * n7 n4 n8
- */
- float p = dog(mln_psite(S)(row, col));
-
- float n1 = dog(mln_psite(S)(row - 1, col));
- float n2 = dog(mln_psite(S)(row, col - 1));
- float n3 = dog(mln_psite(S)(row, col + 1));
- float n4 = dog(mln_psite(S)(row + 1, col));
-
- float n5 = dog(mln_psite(S)(row - 1, col - 1));
- float n6 = dog(mln_psite(S)(row - 1, col + 1));
- float n7 = dog(mln_psite(S)(row + 1, col - 1));
- float n8 = dog(mln_psite(S)(row + 1, col + 1));
-
- // Set Hessian matrix and find the ratio
-
- float dxx = n3 - 2.0 * p + n2;
- float dyy = n4 - 2.0 * p + n1;
- float dxy = (n8 - n6 - n7 + n5) / 4.0;
-
- float trace = dxx + dyy;
- float determinant = dxx * dyy - dxy * dxy;
- float ratio = (trace * trace) / determinant;
-
- if (determinant <= 0 || ratio >= bound_ratio)
- keypoints.erase(k);
+ enough = true;
}
+ }
}
- virtual void save()
+ if (count >= max_interpolation_depth)
{
+ k = keypoints.erase(k);
}
-
- private:
- ScaleSpace<I>* ss;
-
- int isExtrema(S& current, value::int_s16& center,
- typename neighb2d::fwd_niter n)
+ else
{
- bool min = true;
- bool max = true;
-
- for_all(n)
- {
- min = (min && center < current(n));
- max = (max && center > current(n));
- }
-
- if (min)
- return -1;
- else if (max)
- return 1;
- else
- return 0;
+ unsigned octave = k->getOctave();
+ unsigned gradient = k->getGradient();
+
+ S dog = ScaleSpace<S>::pyramid[octave][gradient];
+ S dog_up = ScaleSpace<S>::pyramid[octave][gradient + 1];
+ S dog_down = ScaleSpace<S>::pyramid[octave][gradient - 1];
+ unsigned col = k->getX();
+ unsigned row = k->getY();
+
+ /* Set the neighbors of the current level:
+ * n1
+ * n2 p n3
+ * n4
+ */
+ float p = dog(mln_psite(S)(row, col));
+
+ float n1 = dog(mln_psite(S)(row - 1, col));
+ float n2 = dog(mln_psite(S)(row, col - 1));
+ float n3 = dog(mln_psite(S)(row, col + 1));
+ float n4 = dog(mln_psite(S)(row + 1, col));
+
+ // Set the neighbors of the upper/lower level:
+ float u = dog_up(mln_psite(S)(row, col));
+ float l = dog_down(mln_psite(S)(row, col));
+
+ Matrix firstOrder(3, 1);
+ firstOrder.initialize();
+
+ /************ WARNING ************/
+ /* Maybe missing /255 */
+ /*********************************/
+
+ firstOrder.set(0, 0, n3 - n2);
+ firstOrder.set(1, 0, n4 - n1);
+ firstOrder.set(2, 0, u - l);
+
+ float product = firstOrder.dotProduct(offset);
+ float contrast = p + product / 2.0;
+
+ if (contrast < 0.03)
+ k = keypoints.erase(k);
+ else
+ k++;
}
- void findMin(value::int_s16& min, S& lower, S& upper,
- typename S::piter& p, typename neighb2d::fwd_niter& n)
+ enough = false;
+ count = 0;
+ }
+ }
+
+ void eliminateEdgeResponses(std::list<Keypoint>& keypoints)
+ {
+ const static float bound_ratio = 12.1;
+ std::list<Keypoint>::iterator k = keypoints.begin();
+ std::list<Keypoint>::iterator end = keypoints.end();
+
+ while (k != end)
+ {
+ unsigned octave = k->getOctave();
+ unsigned gradient = k->getGradient();
+ S dog = ScaleSpace<S>::pyramid[octave][gradient];
+ int col = k->getX();
+ int row = k->getY();
+
+ if (row < 0 || col < 0
+ || static_cast<unsigned>(row) >= dog.nrows()
+ || static_cast<unsigned>(col) >= dog.ncols())
{
- min = (lower(p) < min) ? lower(p) : min;
- min = (upper(p) < min) ? upper(p) : min;
-
- for_all(n)
- {
- min = (lower(n) < min) ? lower(n) : min;
- min = (upper(n) < min) ? upper(n) : min;
- }
+ k = keypoints.erase(k);
}
-
- void findMax(value::int_s16& max, S& lower, S& upper,
- typename S::piter& p, typename neighb2d::fwd_niter& n)
+ else
{
- max = (lower(p) > max) ? lower(p) : max;
- max = (upper(p) > max) ? upper(p) : max;
-
- for_all(n)
- {
- max = (lower(n) > max) ? lower(n) : max;
- max = (upper(n) > max) ? upper(n) : max;
- }
+ /* Set the neighbors:
+ * n5 n1 n6
+ * n2 p n3
+ * n7 n4 n8
+ */
+ float p = dog(mln_psite(S)(row, col));
+
+ float n1 = dog(mln_psite(S)(row - 1, col));
+ float n2 = dog(mln_psite(S)(row, col - 1));
+ float n3 = dog(mln_psite(S)(row, col + 1));
+ float n4 = dog(mln_psite(S)(row + 1, col));
+
+ float n5 = dog(mln_psite(S)(row - 1, col - 1));
+ float n6 = dog(mln_psite(S)(row - 1, col + 1));
+ float n7 = dog(mln_psite(S)(row + 1, col - 1));
+ float n8 = dog(mln_psite(S)(row + 1, col + 1));
+
+ // Set Hessian matrix and find the ratio
+
+ float dxx = n3 - 2.0 * p + n2;
+ float dyy = n4 - 2.0 * p + n1;
+ float dxy = (n8 - n6 - n7 + n5) / 4.0;
+
+ float trace = dxx + dyy;
+ float determinant = dxx * dyy - dxy * dxy;
+ float ratio = (trace * trace) / determinant;
+
+ if (determinant <= 0 || ratio >= bound_ratio)
+ k = keypoints.erase(k);
+ else
+ ++k;
}
+ }
+ }
+
+ virtual void save()
+ {
+ ScaleSpace<S>::save("dogs-o");
+ }
+
+ private:
+ ScaleSpace<I>* ss_;
+
+ int isExtrema(S& current, value::int_s16& center,
+ typename neighb2d::fwd_niter n)
+ {
+ bool min = true;
+ bool max = true;
+
+ for_all(n)
+ {
+ min = (min && center < current(n));
+ max = (max && center > current(n));
+ }
+
+ if (min)
+ return -1;
+ else if (max)
+ return 1;
+ else
+ return 0;
+ }
+
+ void findMin(value::int_s16& min, S& lower, S& upper,
+ typename S::piter& p, typename neighb2d::fwd_niter& n)
+ {
+ min = (lower(p) < min) ? lower(p) : min;
+ min = (upper(p) < min) ? upper(p) : min;
+
+ for_all(n)
+ {
+ min = (lower(n) < min) ? lower(n) : min;
+ min = (upper(n) < min) ? upper(n) : min;
+ }
+ }
+
+ void findMax(value::int_s16& max, S& lower, S& upper,
+ typename S::piter& p, typename neighb2d::fwd_niter& n)
+ {
+ max = (lower(p) > max) ? lower(p) : max;
+ max = (upper(p) > max) ? upper(p) : max;
+
+ for_all(n)
+ {
+ max = (lower(n) > max) ? lower(n) : max;
+ max = (upper(n) > max) ? upper(n) : max;
+ }
+ }
};
#endif /* ! DOG_SPACE_HH */
diff --git a/milena/sandbox/anthony/src/main.cc b/milena/sandbox/anthony/src/main.cc
index 5cf4860..a9dd1df 100644
--- a/milena/sandbox/anthony/src/main.cc
+++ b/milena/sandbox/anthony/src/main.cc
@@ -22,58 +22,52 @@ int main(void)
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;
+ std::string source("images/tshirt.pgm");
+ const unsigned blur_level = 3;
+ const unsigned octave_level = 6;
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;
+ C extrema1, extrema1_b, extrema2, extrema2_b, extrema3, extrema3_b;
io::pgm::load(original, source.c_str());
- initialize(extrema_b, original);
+ initialize(extrema1_b, original);
initialize(extrema2_b, original);
- //initialize(improved_b, original);
+ initialize(extrema3_b, original);
- extrema = data::convert(value::rgb8(), original);
+ extrema1 = data::convert(value::rgb8(), original);
extrema2 = data::convert(value::rgb8(), original);
- //improved = data::convert(value::rgb8(), original);
+ extrema3 = data::convert(value::rgb8(), original);
// Localization
ss->build(original);
- //ss->save();
+ ss->save();
dogs->build();
- //dogs->save();
- dogs->findKeypoints(keypoints);
+ dogs->save();
- writeKeypoints(keypoints, extrema, extrema_b);
+ dogs->findKeypoints(keypoints);
+ writeKeypoints(keypoints, extrema1, extrema1_b);
+ std::cout << "Keypoints.size " << keypoints.size() << std::endl;
- std::cout << "BEFORE: Low contrast discarding: " << keypoints.size() << std::endl;
dogs->discardLowContrastKeypoints(keypoints);
- std::cout << "AFTER: Low contrast discarding: " << keypoints.size() << std::endl;
-
- std::cout << "BEFORE: Elimination of edge responses: " << keypoints.size() << std::endl;
- //dogs->eliminateEdgeResponses(keypoints);
- std::cout << "AFTER: Elimination of edge responses: " << keypoints.size() << std::endl;
-
writeKeypoints(keypoints, extrema2, extrema2_b);
+ std::cout << "Keypoints.size " << keypoints.size() << std::endl;
- //writeKeypoints(keypoints, extrema2, extrema2_b);
+ dogs->eliminateEdgeResponses(keypoints);
+ writeKeypoints(keypoints, extrema3, extrema3_b);
+ std::cout << "Keypoints.size " << keypoints.size() << std::endl;
// Save
- io::ppm::save(extrema, "output/extrema.ppm");
- io::ppm::save(extrema_b, "output/extrema_b.ppm");
+ io::ppm::save(extrema1, "output/extrema1.ppm");
+ io::ppm::save(extrema1_b, "output/extrema1_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");
+ io::ppm::save(extrema3, "output/extrema3.ppm");
+ io::ppm::save(extrema3_b, "output/extrema3_b.ppm");
return 0;
}
diff --git a/milena/sandbox/anthony/src/scale-space.hh b/milena/sandbox/anthony/src/scale-space.hh
index e2262b6..d50c07e 100644
--- a/milena/sandbox/anthony/src/scale-space.hh
+++ b/milena/sandbox/anthony/src/scale-space.hh
@@ -10,6 +10,7 @@
#include <mln/io/ppm/save.hh>
#include <mln/core/alias/neighb2d.hh>
#include <mln/literal/all.hh>
+#include <mln/linear/gaussian.hh>
#include <mln/value/int_u8.hh>
#include <mln/value/int_s16.hh>
@@ -62,7 +63,7 @@ class ScaleSpace
}
}
- void save()
+ void save(std::string prefix = "ss-o")
{
std::stringstream name;
@@ -70,7 +71,7 @@ class ScaleSpace
{
for (unsigned j = 0; j < pyramid[i].size(); ++j)
{
- name << "output/ss-o" << i << "-" << j << ".pgm";
+ name << "output/" << prefix << i << "-" << j << ".pgm";
io::pgm::save(pyramid[i][j], name.str());
name.str("");
}
@@ -92,15 +93,19 @@ class ScaleSpace
void addOctave(const I& original)
{
- I blured;
std::vector<I> octave;
- float variance = sqrt(2);
+ float sigma = 1.6;
+ float k = pow(2., 1. / g);
- for (unsigned i = 0; i < g; ++i)
+ octave.push_back(linear::gaussian(original, sigma));
+
+ for (unsigned i = 1; i < g + 3; ++i)
{
- blured = blur(original, variance);
- variance *= sqrt(2);
- octave.push_back(blured);
+ float sig_prev = pow(k, static_cast<float>(i-1)) * sigma;
+ float sig_total = sig_prev * k;
+ float variance = sqrt(sig_total * sig_total - sig_prev * sig_prev);
+
+ octave.push_back(linear::gaussian(original, variance));
}
pyramid.push_back(octave);
--
1.7.10.4
