--- milena/sandbox/anthony/Makefile | 4 +- milena/sandbox/anthony/src/scale_space.cc | 171 ++++++++++++++++++++++++++--- 2 files changed, 157 insertions(+), 18 deletions(-) diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/Makefile index 558bef5..e43db7f 100644 --- a/milena/sandbox/anthony/Makefile +++ b/milena/sandbox/anthony/Makefile @@ -9,10 +9,10 @@ OUTPUT=a.out all: scale -scale: +scale: clean $(CCACHE) $(CC) $(CFLAGS) -O3 -DNDEBUG $(CLIBS) $(SRC) -o $(OUTPUT) -debug: +debug: mrproper $(CCACHE) $(CC) $(CFLAGS) -g -ggdb $(CLIBS) $(SRC) -o $(OUTPUT) clean: diff --git a/milena/sandbox/anthony/src/scale_space.cc b/milena/sandbox/anthony/src/scale_space.cc index 61308c4..0eca2b6 100644 --- a/milena/sandbox/anthony/src/scale_space.cc +++ b/milena/sandbox/anthony/src/scale_space.cc @@ -1,6 +1,9 @@ #include <mln/core/image/image2d.hh> #include <mln/data/all.hh> #include <mln/io/pgm/all.hh> +#include <mln/io/ppm/save.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/literal/all.hh> #include <cmath> @@ -36,13 +39,8 @@ void convolve(const I& original, I& filtered, double *kernel, unsigned t) window2d win; for (int i = -(size / 2); i < (size / 2) + 1; ++i) - { for (int j = -(size / 2); j < (size / 2) + 1; ++j) - { win.insert(i, j); - } - } - // Iterate through all image sites for_all(p) @@ -139,7 +137,7 @@ void buildOctave(const I& original, I blured; initialize(blured, original); - name << "output/o" << n << "b0.pgm"; + name << "output/SS-o" << n << "-0.pgm"; io::pgm::save(original, name.str()); octave.push_back(original); @@ -150,49 +148,190 @@ void buildOctave(const I& original, octave.push_back(blured); name.str(""); - name << "output/o" << n << "b" << i << ".pgm"; + name << "output/SS-o" << n << "-" << i << ".pgm"; io::pgm::save(blured, name.str()); } scaleSpace.push_back(octave); } +// Build n octaves on b blur levels template<typename I> void buildScaleSpace(std::vector< std::vector<I> >& scaleSpace, const I& original, - unsigned octave_level, - unsigned blur_level) + unsigned n, + unsigned b) { - if (octave_level > 2) + if (n > 2) { // 1st octave - buildOctave(original, 1, blur_level, scaleSpace); + buildOctave(original, 1, b, scaleSpace); // 2nd octave I reduced = downscaleResolution(original); - buildOctave(reduced, 2, blur_level, scaleSpace); + buildOctave(reduced, 2, b, scaleSpace); - for (unsigned i = 3; i <= octave_level; ++i) + for (unsigned i = 3; i <= n; ++i) { // i-th octave reduced = downscaleResolution(reduced); - buildOctave(reduced, i, blur_level, scaleSpace); + buildOctave(reduced, i, b, scaleSpace); + } + } + else + { + std::cerr << "ERROR: buildScaleSpace: Wrong level number (" << n << ")" + << std::endl; + } +} + +// Build the DoG pyramid of the scale space +template<typename I> +void buildDifferenceOfGaussianSpace(std::vector< std::vector<I> >& scaleSpace, + std::vector< std::vector<I> >& dogSpace) +{ + for (unsigned i = 0; i < scaleSpace.size(); ++i) + { + std::vector<I> octave = scaleSpace.at(i); + std::vector<I> dogOctave; + + for (unsigned j = 0; j < octave.size() - 1; ++j) + { + std::stringstream name; + I level1 = octave.at(j); + I level2 = octave.at(j + 1); + I dog; + + initialize(dog, level1); + + mln_piter(I) p(level1.domain()); + + for_all(p) + dog(p) = level1(p) - level2(p); + + dogOctave.push_back(dog); + name << "output/DoG_o" << i << "-" << j << ".pgm"; + io::pgm::save(dog, name.str()); } + + dogSpace.push_back(dogOctave); + } +} + +// Find the extremum value in the c8 neighborhood +template<typename T, typename I, typename N> +void find_extremum(T& min, T& max, const I& current, N n) +{ + for_all(n) + { + if (current(n) > max) + max = current(n); + + if (current(n) < min) + min = current(n); } } +// Find the minimum value in the upper and lower layers around p +template<typename T, typename I, typename P> +void find_min(T& min, const I& upper, const I& lower, const P& p) +{ + mln_niter(neighb2d) n(c8(), p); + + for_all(n) + { + if (lower(n) < min) + min = lower(n); + + if (upper(n) < min) + min = upper(n); + } +} + +// Find the maximum value in the upper and lower layers around p +template<typename T, typename I, typename P> +void find_max(T& max, const I& upper, const I& lower, const P& p) +{ + mln_niter(neighb2d) n(c8(), p); + + for_all(n) + { + if (lower(n) > max) + max = lower(n); + + if (upper(n) > max) + max = upper(n); + } +} + +// Build the extrema image +template<typename I, typename C> +void buildExtrema(C extrema, + I original, + std::vector< std::vector<I> >& dogSpace) +{ + extrema = data::convert(value::rgb8(), original); + value::int_u8 min = 255; + value::int_u8 max = 0; + value::int_u8 center = 0; + + for (unsigned i = 0; i < dogSpace.size(); ++i) + { + std::vector<I> octave = dogSpace.at(i); + + for (unsigned j = 1; j < octave.size() - 1; ++j) + { + I current = octave.at(j); + I upper = octave.at(j+1); + I lower = octave.at(j-1); + + mln_piter(I) p(current.domain()); + + for_all(p) + { + mln_niter(neighb2d) n(c8(), p); + center = current(p); + + find_extremum(min, max, current, n); + + if (center <= min) + { + find_min(min, upper, lower, p); + + if (center <= min) + extrema(p) = literal::blue; + } + else if (center >= max) + { + find_max(max, upper, lower, p); + + if (center >= max) + extrema(p) = literal::red; + } + } + } + } + io::ppm::save(extrema, "output/extrema.ppm"); +} + // Main entry point int main(int argc, char** argv) { typedef image2d<value::int_u8> I; - const unsigned blur_level = 3; + typedef image2d<value::rgb8> C; + + const unsigned blur_level = 4; const unsigned octave_level = 4; std::vector< std::vector<I> > scaleSpace; + std::vector< std::vector<I> > dogSpace; I original; + C extrema; - io::pgm::load(original, "images/flower.pgm"); + io::pgm::load(original, "images/keith512.pgm"); buildScaleSpace(scaleSpace, original, octave_level, blur_level); + buildDifferenceOfGaussianSpace(scaleSpace, dogSpace); + buildExtrema(extrema, original, dogSpace); return 0; } -- 1.7.2.5