Working version of gaussian gradient magnitude

src/gaussian_gradient_magnitude.cpp

@@ -2,12 +2,10 @@
 #include <Rcpp.h>
 #include <vigra/multi_array.hxx>
 #include <vigra/multi_convolution.hxx>
-#include <vigra/numerictraits.hxx>
-//#include <vigra/multi_blockwise.hxx>
+#include <vigra/convolution.hxx>
 using namespace Rcpp;
 using namespace vigra;
 
-
 //' Gaussian gradient magnitude
 //'
 //' Magnitude of Gaussian gradient.
@@ -20,11 +18,22 @@ NumericMatrix gaussian_gradient_magnitude(NumericMatrix x, double sigma) {
   const int nrows = x.nrow();
   const int ncols = x.ncol();
   Shape2 shape(nrows, ncols);
-  MultiArrayView<2, double> input_image(shape, x.begin());
+  MultiArrayView<2, double> src(shape, x.begin());
   NumericMatrix result(nrows, ncols);
-  MultiArrayView<2, double> output_image(shape, result.begin());
-  // gaussianGradientMagnitude(input_image, output_image, sigma);
-  //std::copy(output_image.begin(), output_image.end(), result.begin());
+  MultiArrayView<2, double> dest(shape, result.begin());
+  // The following line does not compile, so we calculate things explicitly
+  // gaussianGradientMagnitude(src, dest, sigma);
+  MultiArray<2, TinyVector<float, 2> > grad(src.shape());
+  gaussianGradient(src, grad, sigma);
+  MultiArrayView<2, float, StridedArrayTag> first = grad.bindElementChannel(0);
+  MultiArrayView<2, float, StridedArrayTag> second = grad.bindElementChannel(1);
+  NumericVector first_result(nrows*ncols);
+  NumericVector second_result(nrows*ncols);
+  std::copy(first.begin(), first.end(), first_result.begin());
+  std::copy(second.begin(), second.end(), second_result.begin());
+  NumericVector result_vector(nrows*ncols);
+  result_vector = first_result * first_result + second_result * second_result;
+  Rcpp::algorithm::sqrt(result_vector.begin(), result_vector.end(), result.begin());
   return(result);
 }