[Olena-patches] olena: olena-2.0-597-gb783446 [SIFT] Start new implementation + Build Scale Space

--- milena/sandbox/anthony/{ => new-implem}/Makefile | 18 +- milena/sandbox/anthony/new-implem/src/main.cc | 63 +++++ .../sandbox/anthony/{ => new-implem}/src/matrix.cc | 0 .../sandbox/anthony/{ => new-implem}/src/matrix.hh | 3 +- .../sandbox/anthony/new-implem/src/scale-space.hh | 106 +++++++ milena/sandbox/anthony/new-implem/src/vrac/vrac.hh | 296 ++++++++++++++++++++ 6 files changed, 475 insertions(+), 11 deletions(-) copy milena/sandbox/anthony/{ => new-implem}/Makefile (53%) create mode 100644 milena/sandbox/anthony/new-implem/src/main.cc copy milena/sandbox/anthony/{ => new-implem}/src/matrix.cc (100%) copy milena/sandbox/anthony/{ => new-implem}/src/matrix.hh (99%) create mode 100644 milena/sandbox/anthony/new-implem/src/scale-space.hh create mode 100644 milena/sandbox/anthony/new-implem/src/vrac/vrac.hh diff --git a/milena/sandbox/anthony/Makefile b/milena/sandbox/anthony/new-implem/Makefile similarity index 53% copy from milena/sandbox/anthony/Makefile copy to milena/sandbox/anthony/new-implem/Makefile index 2953013..a4f63d9 100644 --- a/milena/sandbox/anthony/Makefile +++ b/milena/sandbox/anthony/new-implem/Makefile @@ -1,15 +1,15 @@ -CCACHE=ccache +CCACHE= CC=g++ -CFLAGS=-Wall -Werror -CLIBS=-I../../ -CLEAN=*.o output/* log +CFLAGS=-Wall -Wextra +CLIBS=-I../../../ +CLEAN=*.o output/* log debug_* -SRC=src/main.cc src/matrix.cc src/keypoint.cc -OUTPUT=a.out +SRC=src/main.cc src/matrix.cc +OUTPUT=new-implem -all: scale +all: new-implem -scale: clean +new-implem: clean $(CCACHE) $(CC) $(CFLAGS) -O3 -DNDEBUG $(CLIBS) $(SRC) -o $(OUTPUT) debug: mrproper @@ -21,4 +21,4 @@ clean: mrproper: clean rm -f $(OUTPUT) -.PHONY: scale clean mrproper +.PHONY: new-implem clean mrproper diff --git a/milena/sandbox/anthony/new-implem/src/main.cc b/milena/sandbox/anthony/new-implem/src/main.cc new file mode 100644 index 0000000..45a36ad --- /dev/null +++ b/milena/sandbox/anthony/new-implem/src/main.cc @@ -0,0 +1,63 @@ +#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 <mln/value/int_s8.hh> + +//#include <cmath> +#include <string> + +#include "scale-space.hh" +//#include "dog-space.hh" +//#include "keypoint.hh" +//#include "matrix.hh" + +int main(void) +{ + typedef image2d<value::int_u8> I; + typedef image2d<value::rgb8> C; + + // General parameters + std::string source("images/keith.pbm"); + const unsigned blur_level = 5; + const unsigned octave_level = 3; + const bool black = false; + + ScaleSpace<I> *ss = new ScaleSpace<I>(octave_level, blur_level); + //DoGSpace<C>* dog = new DogSpace(); + //std::vector<Keypoint>* keypoints = new std::vector<Keypoint>(); + + I original; + C extrema, improved; + + io::pgm::load(original, source.c_str()); + + if (black) + { + initialize(extrema, original); + initialize(improved, original); + } + else + { + extrema = data::convert(value::rgb8(), original); + improved = data::convert(value::rgb8(), original); + } + + // Localization + ss->build(original); + ss->save(); + //buildDifferenceOfGaussianSpace(scaleSpace, dogSpace); + //buildExtrema(extrema, dogSpace, keypoints); + //discardLowContrastKeypoints(dogSpace, keypoints, improved); + + // Processing + // TODO + + // Save + 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/new-implem/src/matrix.cc similarity index 100% copy from milena/sandbox/anthony/src/matrix.cc copy to milena/sandbox/anthony/new-implem/src/matrix.cc diff --git a/milena/sandbox/anthony/src/matrix.hh b/milena/sandbox/anthony/new-implem/src/matrix.hh similarity index 99% copy from milena/sandbox/anthony/src/matrix.hh copy to milena/sandbox/anthony/new-implem/src/matrix.hh index de83912..c42de59 100644 --- a/milena/sandbox/anthony/src/matrix.hh +++ b/milena/sandbox/anthony/new-implem/src/matrix.hh @@ -16,6 +16,7 @@ class Matrix bool inverse(Matrix& inversed); void absolute(); bool isMajorOffset(); + float determinant(); // Operators Matrix operator*(Matrix& right); @@ -31,8 +32,6 @@ class Matrix std::vector< std::vector<float> > m; unsigned height; unsigned width; - - float determinant(); }; #endif /* ! MATRIX_HH */ diff --git a/milena/sandbox/anthony/new-implem/src/scale-space.hh b/milena/sandbox/anthony/new-implem/src/scale-space.hh new file mode 100644 index 0000000..2fbb4fd --- /dev/null +++ b/milena/sandbox/anthony/new-implem/src/scale-space.hh @@ -0,0 +1,106 @@ +#ifndef SCALE_SPACE_HH +# define SCALE_SPACE_HH + +# include <string> +# include <vector> + +#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 <mln/value/int_s8.hh> + +#include "vrac/vrac.hh" + +using namespace mln; + +template<typename I> +class ScaleSpace +{ + public: + ScaleSpace(unsigned octave_level, unsigned gradient_level) + { + o = octave_level; + g = gradient_level; + } + + ~ScaleSpace() + { + typedef typename std::vector< std::vector<I> >::iterator PIT; + + for (PIT it = pyramid.begin(); it != pyramid.end(); ++it) + it->clear(); + + pyramid.clear(); + } + + std::vector<I>& getOctave(unsigned i) + { + return pyramid[i]; + } + + + void build(const I& original) + { + // Upscaled octave + I upscaled = vrac::upscaling2(original); + addOctave(upscaled); + + // 1st octave + addOctave(original); + + // 2nd octave + I downscaled = vrac::downscaling2(original); + addOctave(downscaled); + + for (unsigned i = 3; i < o; ++i) + { + downscaled = vrac::downscaling2(downscaled); + addOctave(downscaled); + } + } + + void save() + { + std::stringstream name; + + for (unsigned i = 0; i < pyramid.size(); ++i) + { + for (unsigned j = 0; j < pyramid[i].size(); ++j) + { + name << "output/ss-o" << i << "-" << j << ".pgm"; + io::pgm::save(pyramid[i][j], name.str()); + name.str(""); + } + } + } + + inline unsigned octave_count() { return pyramid.size(); } + inline unsigned gradient_count() { return pyramid[0].size(); } + + private: + unsigned o; + unsigned g; + std::vector< std::vector<I> > pyramid; + + void addOctave(const I& original) + { + I blured; + std::vector<I> octave; + //float variance = sqrt(2); + float variance = 0.84089642; + + for (unsigned i = 0; i < g; ++i) + { + blured = vrac::blur(original, variance); + variance *= sqrt(2); + octave.push_back(blured); + } + + pyramid.push_back(octave); + } +}; + +#endif /* ! SCALE_SPACE_HH */ diff --git a/milena/sandbox/anthony/new-implem/src/vrac/vrac.hh b/milena/sandbox/anthony/new-implem/src/vrac/vrac.hh new file mode 100644 index 0000000..fc9e3f2 --- /dev/null +++ b/milena/sandbox/anthony/new-implem/src/vrac/vrac.hh @@ -0,0 +1,296 @@ +#ifndef VRAC_HH +# define VRAC_HH + +# 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 <mln/value/int_s8.hh> + +# include <cmath> + +# include "../matrix.hh" + +//# include "keypoint.hh" +//# include "matrix.hh" + +using namespace mln; + +namespace vrac +{ + // Set the gaussian kernel + float *gaussian_kernel(float variance) + { + unsigned index = 0; + float *kernel = new float[49]; + Matrix matrix(7, 7); + float sum = 0; + + for (int i = -3; i < 4; ++i) + { + for (int j = -3; j < 4; ++j) + { + kernel[index] = exp(- (i*i + j*j) / (2.0 * variance * variance) ); + sum += kernel[index++]; + } + } + + sum = 1. / sum; + + for (int i = 0; i < 49; ++i) + kernel[i] *= sum; + + return kernel; + } + + // Apply the convolution of the image by the kernel + template<typename I> + void convolve(const I& original, I& filtered, float *kernel) + { + mln_piter(I) p(filtered.domain()); + window2d win; + + for (int i = -3; i < 4; ++i) + for (int j = -3; j < 4; ++j) + win.insert(i, j); + + // Iterate through all image sites + for_all(p) + { + // Create the window around the site + mln_qiter(window2d) q(win, p); + float sum = 0; + int index = 0; + + // Iterate through all window image sites + for_all(q) + { + if (filtered.has(q)) + sum += original(q) * kernel[index++]; + else + sum += 127 * kernel[index++]; + } + + filtered(p) = static_cast<int>(sum); + } + } + + // Blur the image thanks to a convolution matrix + template<typename I> + I blur(const I& original, float variance) + { + std::cout << "variance: " << variance << std::endl; + I filtered; + float *kernel = gaussian_kernel(variance); + + initialize(filtered, original); + convolve(original, filtered, kernel); + + return filtered; + } + + float *upscalingTL_kernel() + { + const static float one = 1.0 / 8.0; + const static float two = 2.0 / 8.0; + const static float three = 2.0 / 8.0; + const static float four = 3.0 / 8.0; + + float *kernel = new float[4]; + + kernel[0] = one; + kernel[1] = two; + kernel[2] = three; + kernel[3] = four; + + return kernel; + } + + float *upscalingTR_kernel() + { + const static float one = 2.0 / 8.0; + const static float two = 1.0 / 8.0; + const static float three = 3.0 / 8.0; + const static float four = 2.0 / 8.0; + + float *kernel = new float[4]; + + kernel[0] = one; + kernel[1] = two; + kernel[2] = three; + kernel[3] = four; + + return kernel; + } + + float *upscalingBL_kernel() + { + const static float one = 2.0 / 8.0; + const static float two = 3.0 / 8.0; + const static float three = 1.0 / 8.0; + const static float four = 2.0 / 8.0; + + float *kernel = new float[4]; + + kernel[0] = one; + kernel[1] = two; + kernel[2] = three; + kernel[3] = four; + + return kernel; + } + + float *upscalingBR_kernel() + { + const static float one = 3.0 / 8.0; + const static float two = 2.0 / 8.0; + const static float three = 2.0 / 8.0; + const static float four = 1.0 / 8.0; + + float *kernel = new float[4]; + + kernel[0] = one; + kernel[1] = two; + kernel[2] = three; + kernel[3] = four; + + return kernel; + } + + template<typename I> + unsigned upscalePixel(const I& original, int i, int j, int type) + { + float *kernel; + window2d win; + mln_psite(I) p(i, j); + float value = 0; + unsigned index = 0; + int iStart = 0; + int jStart = 0; + + if (type == 0) + { + kernel = upscalingTL_kernel(); + iStart = -1; + jStart = -1; + } + else if (type == 1) + { + kernel = upscalingTR_kernel(); + iStart = 0; + jStart = -1; + } + else if (type == 2) + { + kernel = upscalingBL_kernel(); + iStart = -1; + jStart = 0; + } + else + { + kernel = upscalingBR_kernel(); + iStart = 0; + jStart = 0; + } + + for (int ii = iStart; ii < iStart + 2; ++ii) + for (int jj = jStart; jj < jStart + 2; ++jj) + win.insert(ii, jj); + + mln_qiter(window2d) q(win, p); + + for_all(q) + { + if (original.has(q)) + { + value += original(q) * kernel[index]; + } + + ++index; + } + + return (static_cast<int>(value)); + } + // Downscale by 2 the image resolution + template<typename I> + I downscaling2(const I& original) + { + mln_piter(I) p(original.domain()); + window2d win; + int index = 0; + + // Initialize the rescaled image + box2d size(original.nrows() / 2, original.ncols() / 2); + I reduced(size); + + // Set the 2x2 window matrix + for (int i = 0; i < 2; ++i) + for (int j = 0; j < 2; ++j) + win.insert(i, j); + + // Iterate through all image sites + for_all(p) + { + // Get 1/4 of total pixels + if (index % 2 == 0) + { + int i = (index % original.ncols()) / 2; + int j = (index / original.ncols()) / 2; + int count = 0; + int sum = 0; + mln_qiter(window2d) q(win, p); + + // Iterate through all window image sites + for_all(q) + { + if (original.has(q)) + { + sum += original(q); + ++count; + } + } + + if (reduced.has(point2d(j, i))) + opt::at(reduced, j, i) = sum / count; + } + + ++index; + } + + return reduced; + } + + // Upscale by 2 the image resolution + template<typename I> + I upscaling2(const I& original) + { + box2d size(original.nrows() * 2, original.ncols() * 2); + I upscaled(size); + + mln_piter(I) p(upscaled.domain()); + int index = 0; + + for_all(p) + { + int i = p.col(); + int j = p.row(); + + if (i % 2 == 0 && j % 2 == 0) // Top left + opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 0); + else if (i % 2 == 0) // Bottom left + opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 2); + else if (j % 2 == 0) // Top right + opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 1); + else // Bottom right + opt::at(upscaled, i, j) = upscalePixel(original, i / 2, j / 2, 3); + + ++index; + } + + return upscaled; + } +} + +#endif /* ! VRAC_HH */ -- 1.7.10.4