[Olena-patches] olena: olena-2.0-604-g586804e [SIFT] Minor fixes

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