Newer
Older
//[[Rcpp::depends(vigRa)]]
#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());
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());