gaussian_gradient_magnitude.cpp

//[[Rcpp::depends(VH)]]
#include <Rcpp.h>
#include <vigra/multi_array.hxx>
#include <vigra/multi_convolution.hxx>
#include <vigra/convolution.hxx>
using namespace Rcpp;
using namespace vigra;

//' Gaussian gradient magnitude
//'
//' Magnitude of Gaussian gradient.
//' @param x Numeric matrix
//' @param sigma Scale of the Gaussian kernel
//' @return Numeric matrix; smoothed version of x
//' @export
// [[Rcpp::export]]
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> src(shape, x.begin());
  NumericMatrix result(nrows, ncols);
  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);
}