--- milena/sandbox/anthony/src/main.cc | 583 ++++++++++++++++++++++++++++++++++++ 1 files changed, 583 insertions(+), 0 deletions(-) create mode 100644 milena/sandbox/anthony/src/main.cc diff --git a/milena/sandbox/anthony/src/main.cc b/milena/sandbox/anthony/src/main.cc new file mode 100644 index 0000000..eea1c86 --- /dev/null +++ b/milena/sandbox/anthony/src/main.cc @@ -0,0 +1,583 @@ +#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 "keypoint.hh" + +using namespace mln; + +// Set the gaussian kernel +float *gaussian_kernel() +{ + float variance = sqrt(2.0); + unsigned index = 0; + float *kernel = new float[49]; + static const float pi = 4.0 * atan(1.0); + float div = (2.0 * pi * variance); + + for (int i = -3; i < 4; ++i) + { + for (int j = -3; j < 4; ++j) + { + float e = exp(- (i*i + j*j) / (2.0 * variance) ); + kernel[index++] = e / div; + } + } + + 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) +{ + I filtered; + float *kernel = gaussian_kernel(); + + initialize(filtered, original); + convolve(original, filtered, kernel); + + return filtered; +} + +// 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; +} + +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)); +} + +// 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; +} + +// Build the n-th octave on b blur levels +template<typename I> +void buildOctave(const I& original, + unsigned n, + unsigned b, + std::vector< std::vector<I> >& scaleSpace) +{ + std::vector<I> octave; + std::stringstream name; + I blured; + initialize(blured, original); + + name << "output/SS-o" << n << "-0.pgm"; + io::pgm::save(original, name.str()); + + octave.push_back(original); + blured = duplicate(original); + + for(unsigned i = 1; i <= b; i *= 2) + { + blured = blur(blured); + octave.push_back(blured); + + name.str(""); + 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 n, + unsigned b) +{ + if (n > 2) + { + // Upscaled octave + I upscaled = upscaling2(original); + buildOctave(upscaled, 0, b, scaleSpace); + + // 1st octave + buildOctave(original, 1, b, scaleSpace); + + // 2nd octave + I downscaled = downscaling2(original); + buildOctave(downscaled, 2, b, scaleSpace); + + for (unsigned i = 3; i <= n; ++i) + { + // i-th octave + downscaled = downscaling2(downscaled); + buildOctave(downscaled, 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) + { + int diff = ((int) level1(p)) - ((int) level2(p)); + dog(p) = (diff >= 0) ? diff : 0; + } + + 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& pUpper, const P& pLower) +{ + mln_niter(neighb2d) nUpper(c8(), pUpper); + mln_niter(neighb2d) nLower(c8(), pLower); + + // Upper processing + for_all(nUpper) + if (upper(nUpper) < min) + min = upper(nUpper); + + // Lower processing + for_all(nLower) + if (lower(nLower) < min) + min = lower(nLower); +} + +// 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& pUpper, const P& pLower) +{ + mln_niter(neighb2d) nUpper(c8(), pUpper); + mln_niter(neighb2d) nLower(c8(), pLower); + + // Upper processing + for_all(nUpper) + if (upper(nUpper) > max) + max = upper(nUpper); + + // Lower processing + for_all(nLower) + if (lower(nLower) > max) + max = lower(nLower); +} + +template<typename T, typename I, typename P> +bool hasEdgeResponse(const I& dog, const P& p, const T& center) +{ + const static float bound_ratio = 12.1; + + /* Set the neighbors: + * n5 n1 n6 + * n2 p n3 + * n7 n4 n8 + */ + float n1 = dog(mln_psite(I)(p.row() - 1, p.col())); + float n2 = dog(mln_psite(I)(p.row(), p.col() - 1)); + float n3 = dog(mln_psite(I)(p.row(), p.col() + 1)); + float n4 = dog(mln_psite(I)(p.row() + 1, p.col())); + + float n5 = dog(mln_psite(I)(p.row() - 1, p.col() - 1)); + float n6 = dog(mln_psite(I)(p.row() - 1, p.col() + 1)); + float n7 = dog(mln_psite(I)(p.row() + 1, p.col() - 1)); + float n8 = dog(mln_psite(I)(p.row() + 1, p.col() + 1)); + + // Set Hessian matrix and find the ratio + + float dxx = n3 - 2.0 * center + n2; + float dyy = n4 - 2.0 * center + n1; + + float dxy = (n8 - n6 - n7 + n5) / 4.0; + float trace = dxx + dyy; + float determinant = dxx * dyy - dxy * dxy; + + if (determinant <= 0) + return true; + + float ratio = (trace * trace) / determinant; + + return (ratio >= bound_ratio); +} + +// Build the extrema image +template<typename I, typename C> +void buildExtrema(C extrema, + 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(); + + // Iterates through all scales + for (unsigned i = 0; i < scales; ++i) + { + // Iterates through all resolutions + for (unsigned j = 1; j < dogSpace.size() - 1; ++j) + { + I current = dogSpace.at(j).at(i); + I upper = dogSpace.at(j+1).at(i); + I lower = dogSpace.at(j-1).at(i); + + mln_piter(I) p(current.domain()); + + for_all(p) + { + value::int_u8 min = 255; + value::int_u8 max = 0; + value::int_u8 center = 0; + + mln_niter(neighb2d) n(c8(), p); + center = current(p); + + mln_psite(I) pUpper(p.row() / 2, p.col() / 2); + mln_psite(I) pLower(p.row() * 2, p.col() * 2); + + // Find the min and max in the current neighborhood (c8) + find_extremum(min, max, current, n); + + if (center <= min) + { + // Find the min in other levels of the pyramid + find_min(min, upper, lower, pUpper, pLower); + + if (center < min && !hasEdgeResponse(current, p, center)) + { + mln_psite(I) pOriginal(p.row() * pow(2, (j-1)), + p.col() * pow(2, (j-1))); + extrema(pOriginal) = literal::green; + keypoints.push_back(Keypoint(p.row(), p.col(), j, i, false)); + } + } + else if (center >= max) + { + // Find the max in other levels of the pyramid + find_max(max, upper, lower, pUpper, pLower); + + if (center > max && !hasEdgeResponse(current, p, center)) + { + mln_psite(I) pOriginal(p.row() * pow(2, (j-1)), + p.col() * pow(2, (j-1))); + extrema(pOriginal) = literal::red; + keypoints.push_back(Keypoint(p.row(), p.col(), j, i, true)); + } + } + } + } + } + + io::ppm::save(extrema, "output/extrema.ppm"); +} + +template<typename I> +void discardEdgeResponses(std::vector< std::vector<I> >& dogSpace, + std::vector<Keypoint>& keypoints) +{ + // FIXME Move work from buildExtrema here +} + +// Main entry point +int main(int argc, char** argv) +{ + typedef image2d<value::int_u8> I; + typedef image2d<value::rgb8> C; + + const unsigned blur_level = 5; + const unsigned octave_level = 3; + std::vector< std::vector<I> > scaleSpace; + std::vector< std::vector<I> > dogSpace; + std::vector<Keypoint> keypoints; + I original; + C extrema; + + io::pgm::load(original, "images/lena.pgm"); + + buildScaleSpace(scaleSpace, original, octave_level, blur_level); + buildDifferenceOfGaussianSpace(scaleSpace, dogSpace); + buildExtrema(extrema, original, dogSpace, keypoints); + // TODO Interpolation of nearby data for accurate position + // TODO Discard low-contrast keypoints + //discardEdgeResponses(dogSpace, keypoints); + + std::cout << keypoints.size() << std::endl; + + + return 0; +} -- 1.7.2.5