---
milena/sandbox/anthony/src/scale_space.cc | 191 ++++++++++++++++++++++++++---
1 files changed, 174 insertions(+), 17 deletions(-)
diff --git a/milena/sandbox/anthony/src/scale_space.cc
b/milena/sandbox/anthony/src/scale_space.cc
index 735d3f1..87853ef 100644
--- a/milena/sandbox/anthony/src/scale_space.cc
+++ b/milena/sandbox/anthony/src/scale_space.cc
@@ -11,19 +11,19 @@
using namespace mln;
// Set the gaussian kernel
-double *gaussian_kernel()
+float *gaussian_kernel()
{
- double variance = sqrt(2.0);
+ float variance = sqrt(2.0);
unsigned index = 0;
- double *kernel = new double[49];
- static const double pi = 4.0 * atan(1.0);
- double div = (2.0 * pi * variance);
+ 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)
{
- double e = exp(- (i*i + j*j) / (2.0 * variance) );
+ float e = exp(- (i*i + j*j) / (2.0 * variance) );
kernel[index++] = e / div;
}
}
@@ -33,7 +33,7 @@ double *gaussian_kernel()
// Apply the convolution of the image by the kernel
template<typename I>
-void convolve(const I& original, I& filtered, double *kernel)
+void convolve(const I& original, I& filtered, float *kernel)
{
mln_piter(I) p(filtered.domain());
window2d win;
@@ -68,7 +68,7 @@ template<typename I>
I blur(const I& original)
{
I filtered;
- double *kernel = gaussian_kernel();
+ float *kernel = gaussian_kernel();
initialize(filtered, original);
convolve(original, filtered, kernel);
@@ -125,6 +125,160 @@ I downscaling2(const I& original)
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,
@@ -165,18 +319,22 @@ void buildScaleSpace(std::vector< std::vector<I> >&
scaleSpace,
{
if (n > 2)
{
+ // Upscaled octave
+ I upscaled = upscaling2(original);
+ buildOctave(upscaled, 0, b, scaleSpace);
+
// 1st octave
buildOctave(original, 1, b, scaleSpace);
// 2nd octave
- I reduced = downscaling2(original);
- buildOctave(reduced, 2, b, scaleSpace);
+ I downscaled = downscaling2(original);
+ buildOctave(downscaled, 2, b, scaleSpace);
for (unsigned i = 3; i <= n; ++i)
{
// i-th octave
- reduced = downscaling2(reduced);
- buildOctave(reduced, i, b, scaleSpace);
+ downscaled = downscaling2(downscaled);
+ buildOctave(downscaled, i, b, scaleSpace);
}
}
else
@@ -279,9 +437,6 @@ void buildExtrema(C extrema,
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;
const unsigned blur_level = dogSpace.at(0).size();
@@ -316,7 +471,8 @@ void buildExtrema(C extrema,
if (center < min)
{
- mln_psite(I) pOriginal(p.row() * pow(2, j), p.col() * pow(2, j));
+ mln_psite(I) pOriginal(p.row() * pow(2, (j-1)),
+ p.col() * pow(2, (j-1)));
extrema(pOriginal) = literal::red;
}
}
@@ -326,7 +482,8 @@ void buildExtrema(C extrema,
if (center > max)
{
- mln_psite(I) pOriginal(p.row() * pow(2, j), p.col() * pow(2, j));
+ mln_psite(I) pOriginal(p.row() * pow(2, (j-1)),
+ p.col() * pow(2, (j-1)));
extrema(pOriginal) = literal::red;
}
}
--
1.7.2.5