diff --git a/NAMESPACE b/NAMESPACE
index 2485ee900e40c5500b78de4d2c885229a74e9bad..d3182ab45cd811359afbfe29d1b96928f98d53c6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+export(difference_of_gaussians)
+export(gaussian_gradient_magnitude)
 export(gaussian_smoothing)
+export(hessian_of_gaussian_eigenvalues)
+export(laplacian_of_gaussian)
+export(structure_tensor_eigenvalues)
+export(structure_tensor_eigenvalues_2)
 importFrom(Rcpp,sourceCpp)
 useDynLib(vigratools, .registration = TRUE)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 5ffb9599c90a6e0b4a2993c1c6f257f7fb902ef7..6e5c77724eccfe8a1097ea27f265c3388817e0f8 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,6 +1,17 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' 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
+gaussian_gradient_magnitude <- function(x, sigma) {
+    .Call(`_vigratools_gaussian_gradient_magnitude`, x, sigma)
+}
+
 #' Gaussian smoothing
 #'
 #' Perform isotropic Gaussian convolution. The function uses BORDER_TREATMENT_REFLECT.
@@ -12,3 +23,48 @@ gaussian_smoothing <- function(x, sigma) {
     .Call(`_vigratools_gaussian_smoothing`, x, sigma)
 }
 
+#' Hessian of Gaussian eigenvalues
+#'
+#' Calculate eigenvalues of Hessian of Gaussian
+#' @param x Numeric matrix
+#' @param sigma Scale of smoothing
+#' @return List of 2 matrices with 1st and 2nd eigenvalue
+#' @export
+hessian_of_gaussian_eigenvalues <- function(x, sigma) {
+    .Call(`_vigratools_hessian_of_gaussian_eigenvalues`, x, sigma)
+}
+
+#' Laplacian of Gaussian
+#'
+#' Calculates laplacian of Gaussian.
+#' @param x Numeric matrix
+#' @param sigma Scale of the Gaussian kernel
+#' @return Numeric matrix; smoothed version of x
+#' @export
+laplacian_of_gaussian <- function(x, sigma) {
+    .Call(`_vigratools_laplacian_of_gaussian`, x, sigma)
+}
+
+#' Structure tensor eigenvalues (2 argument)
+#'
+#' Calculate eigenvalues of structure tensor
+#' @param x Numeric matrix
+#' @param inner_sigma Scale of inner smoothing
+#' @param outer_sigma Scale of outer smoothing
+#' @return List of 2 matrices with 1st and 2nd eigenvalue
+#' @export
+structure_tensor_eigenvalues_2 <- function(x, inner_sigma, outer_sigma) {
+    .Call(`_vigratools_structure_tensor_eigenvalues_2`, x, inner_sigma, outer_sigma)
+}
+
+#' Structure tensor eigenvalues
+#'
+#' Calculate eigenvalues of structure tensor
+#' @param x Numeric matrix
+#' @param sigma Scale of inner smoothing; outer smoothing is sigma/2.
+#' @return List of 2 matrices with 1st and 2nd eigenvalue
+#' @export
+structure_tensor_eigenvalues <- function(x, sigma) {
+    .Call(`_vigratools_structure_tensor_eigenvalues`, x, sigma)
+}
+
diff --git a/R/difference_of_gaussians.R b/R/difference_of_gaussians.R
new file mode 100644
index 0000000000000000000000000000000000000000..ee7746ab63717374760916dd9eb80bb338b8c033
--- /dev/null
+++ b/R/difference_of_gaussians.R
@@ -0,0 +1,13 @@
+#' Difference of Gaussians
+#'
+#' Compute difference of Gaussians
+#' @param x Input matrix
+#' @param sigma Scale of smoothing (smaller gaussian is 0.66 times sigma)
+#' @return Matrix containing difference of original image smoothed at two scales.
+#' @export
+#'
+difference_of_gaussians <- function(x, sigma) {
+  smoothed_1 <- gaussian_smoothing(x, sigma)
+  smoothed_2 <- gaussian_smoothing(x, 0.66*sigma)
+  smoothed_1 - smoothed_2
+}
diff --git a/man/difference_of_gaussians.Rd b/man/difference_of_gaussians.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..a1066fb84decbd27785f5f498b31e7821cd194e6
--- /dev/null
+++ b/man/difference_of_gaussians.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/difference_of_gaussians.R
+\name{difference_of_gaussians}
+\alias{difference_of_gaussians}
+\title{Difference of Gaussians}
+\usage{
+difference_of_gaussians(x, sigma)
+}
+\arguments{
+\item{x}{Input matrix}
+
+\item{sigma}{Scale of smoothing (smaller gaussian is 0.66 times sigma)}
+}
+\value{
+Matrix containing difference of original image smoothed at two scales.
+}
+\description{
+Compute difference of Gaussians
+}
diff --git a/man/gaussian_gradient_magnitude.Rd b/man/gaussian_gradient_magnitude.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..5c22da45ac66b36a8642d3561fcd801249570da0
--- /dev/null
+++ b/man/gaussian_gradient_magnitude.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{gaussian_gradient_magnitude}
+\alias{gaussian_gradient_magnitude}
+\title{Gaussian gradient magnitude}
+\usage{
+gaussian_gradient_magnitude(x, sigma)
+}
+\arguments{
+\item{x}{Numeric matrix}
+
+\item{sigma}{Scale of the Gaussian kernel}
+}
+\value{
+Numeric matrix; smoothed version of x
+}
+\description{
+Magnitude of Gaussian gradient.
+}
diff --git a/man/hessian_of_gaussian_eigenvalues.Rd b/man/hessian_of_gaussian_eigenvalues.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..6c7d075bcf1fa6e88db5fbc54267b3ea1143b256
--- /dev/null
+++ b/man/hessian_of_gaussian_eigenvalues.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{hessian_of_gaussian_eigenvalues}
+\alias{hessian_of_gaussian_eigenvalues}
+\title{Hessian of Gaussian eigenvalues}
+\usage{
+hessian_of_gaussian_eigenvalues(x, sigma)
+}
+\arguments{
+\item{x}{Numeric matrix}
+
+\item{sigma}{Scale of smoothing}
+}
+\value{
+List of 2 matrices with 1st and 2nd eigenvalue
+}
+\description{
+Calculate eigenvalues of Hessian of Gaussian
+}
diff --git a/man/laplacian_of_gaussian.Rd b/man/laplacian_of_gaussian.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..3a7c47d4f322f428d7ecdf1fa6475f8dd6555aca
--- /dev/null
+++ b/man/laplacian_of_gaussian.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{laplacian_of_gaussian}
+\alias{laplacian_of_gaussian}
+\title{Laplacian of Gaussian}
+\usage{
+laplacian_of_gaussian(x, sigma)
+}
+\arguments{
+\item{x}{Numeric matrix}
+
+\item{sigma}{Scale of the Gaussian kernel}
+}
+\value{
+Numeric matrix; smoothed version of x
+}
+\description{
+Calculates laplacian of Gaussian.
+}
diff --git a/man/structure_tensor_eigenvalues.Rd b/man/structure_tensor_eigenvalues.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..5f43fdab893a4938b806cc9e9148cab09d2ff409
--- /dev/null
+++ b/man/structure_tensor_eigenvalues.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{structure_tensor_eigenvalues}
+\alias{structure_tensor_eigenvalues}
+\title{Structure tensor eigenvalues}
+\usage{
+structure_tensor_eigenvalues(x, sigma)
+}
+\arguments{
+\item{x}{Numeric matrix}
+
+\item{sigma}{Scale of inner smoothing; outer smoothing is sigma/2.}
+}
+\value{
+List of 2 matrices with 1st and 2nd eigenvalue
+}
+\description{
+Calculate eigenvalues of structure tensor
+}
diff --git a/man/structure_tensor_eigenvalues_2.Rd b/man/structure_tensor_eigenvalues_2.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..94712cc41f3928cd27a4bbce104d730526d42d49
--- /dev/null
+++ b/man/structure_tensor_eigenvalues_2.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{structure_tensor_eigenvalues_2}
+\alias{structure_tensor_eigenvalues_2}
+\title{Structure tensor eigenvalues (2 argument)}
+\usage{
+structure_tensor_eigenvalues_2(x, inner_sigma, outer_sigma)
+}
+\arguments{
+\item{x}{Numeric matrix}
+
+\item{inner_sigma}{Scale of inner smoothing}
+
+\item{outer_sigma}{Scale of outer smoothing}
+}
+\value{
+List of 2 matrices with 1st and 2nd eigenvalue
+}
+\description{
+Calculate eigenvalues of structure tensor
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 9301353b1bf6069c722a0c9494802eb60d769dcd..0d71229e27f7dfef07fc685f7024d851eb72a446 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -5,6 +5,18 @@
 
 using namespace Rcpp;
 
+// gaussian_gradient_magnitude
+NumericMatrix gaussian_gradient_magnitude(NumericMatrix x, double sigma);
+RcppExport SEXP _vigratools_gaussian_gradient_magnitude(SEXP xSEXP, SEXP sigmaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
+    rcpp_result_gen = Rcpp::wrap(gaussian_gradient_magnitude(x, sigma));
+    return rcpp_result_gen;
+END_RCPP
+}
 // gaussian_smoothing
 NumericMatrix gaussian_smoothing(NumericMatrix x, double sigma);
 RcppExport SEXP _vigratools_gaussian_smoothing(SEXP xSEXP, SEXP sigmaSEXP) {
@@ -17,9 +29,63 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// hessian_of_gaussian_eigenvalues
+List hessian_of_gaussian_eigenvalues(NumericMatrix x, double sigma);
+RcppExport SEXP _vigratools_hessian_of_gaussian_eigenvalues(SEXP xSEXP, SEXP sigmaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
+    rcpp_result_gen = Rcpp::wrap(hessian_of_gaussian_eigenvalues(x, sigma));
+    return rcpp_result_gen;
+END_RCPP
+}
+// laplacian_of_gaussian
+NumericMatrix laplacian_of_gaussian(NumericMatrix x, double sigma);
+RcppExport SEXP _vigratools_laplacian_of_gaussian(SEXP xSEXP, SEXP sigmaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
+    rcpp_result_gen = Rcpp::wrap(laplacian_of_gaussian(x, sigma));
+    return rcpp_result_gen;
+END_RCPP
+}
+// structure_tensor_eigenvalues_2
+List structure_tensor_eigenvalues_2(NumericMatrix x, double inner_sigma, double outer_sigma);
+RcppExport SEXP _vigratools_structure_tensor_eigenvalues_2(SEXP xSEXP, SEXP inner_sigmaSEXP, SEXP outer_sigmaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type inner_sigma(inner_sigmaSEXP);
+    Rcpp::traits::input_parameter< double >::type outer_sigma(outer_sigmaSEXP);
+    rcpp_result_gen = Rcpp::wrap(structure_tensor_eigenvalues_2(x, inner_sigma, outer_sigma));
+    return rcpp_result_gen;
+END_RCPP
+}
+// structure_tensor_eigenvalues
+List structure_tensor_eigenvalues(NumericMatrix x, double sigma);
+RcppExport SEXP _vigratools_structure_tensor_eigenvalues(SEXP xSEXP, SEXP sigmaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
+    rcpp_result_gen = Rcpp::wrap(structure_tensor_eigenvalues(x, sigma));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_vigratools_gaussian_gradient_magnitude", (DL_FUNC) &_vigratools_gaussian_gradient_magnitude, 2},
     {"_vigratools_gaussian_smoothing", (DL_FUNC) &_vigratools_gaussian_smoothing, 2},
+    {"_vigratools_hessian_of_gaussian_eigenvalues", (DL_FUNC) &_vigratools_hessian_of_gaussian_eigenvalues, 2},
+    {"_vigratools_laplacian_of_gaussian", (DL_FUNC) &_vigratools_laplacian_of_gaussian, 2},
+    {"_vigratools_structure_tensor_eigenvalues_2", (DL_FUNC) &_vigratools_structure_tensor_eigenvalues_2, 3},
+    {"_vigratools_structure_tensor_eigenvalues", (DL_FUNC) &_vigratools_structure_tensor_eigenvalues, 2},
     {NULL, NULL, 0}
 };
 
diff --git a/src/gaussian_gradient_magnitude.cpp b/src/gaussian_gradient_magnitude.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..73f9a190a01cec80c171779259d81730ef5aa801
--- /dev/null
+++ b/src/gaussian_gradient_magnitude.cpp
@@ -0,0 +1,31 @@
+//[[Rcpp::depends(vigRa)]]
+#include <Rcpp.h>
+#include <vigra/multi_array.hxx>
+#include <vigra/multi_convolution.hxx>
+#include <vigra/numerictraits.hxx>
+//#include <vigra/multi_blockwise.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> input_image(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());
+  return(result);
+}
+
+
diff --git a/src/gaussian_smoothing.cpp b/src/gaussian_smoothing.cpp
index df3ea6cb20e611dd92dfae5a1b5fd55dbf38bdca..a1fe864f35fd2866a9a3828b2bc7a471a9b2ad38 100644
--- a/src/gaussian_smoothing.cpp
+++ b/src/gaussian_smoothing.cpp
@@ -2,9 +2,6 @@
 #include <Rcpp.h>
 #include <vigra/multi_array.hxx>
 #include <vigra/convolution.hxx>
-// #include <vigra/tensorutilities.hxx>
-// #include <vigra/mathutil.hxx>
-// #include <vigra/multi_tensorutilities.hxx>
 using namespace Rcpp;
 using namespace vigra;
 
@@ -27,3 +24,5 @@ NumericMatrix gaussian_smoothing(NumericMatrix x, double sigma) {
   gaussianSmoothing(input_image, output_image, sigma);
   return(result);
 }
+
+
diff --git a/src/hessian_of_gaussian_eigenvalues.cpp b/src/hessian_of_gaussian_eigenvalues.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f947b516a49640c79656ab9d478b7c33e8466fae
--- /dev/null
+++ b/src/hessian_of_gaussian_eigenvalues.cpp
@@ -0,0 +1,36 @@
+//[[Rcpp::depends(vigRa)]]
+#include <Rcpp.h>
+#include <vigra/multi_array.hxx>
+#include <vigra/convolution.hxx>
+#include <vigra/tensorutilities.hxx>
+#include <vigra/mathutil.hxx>
+#include <vigra/multi_tensorutilities.hxx>
+using namespace Rcpp;
+using namespace vigra;
+
+//' Hessian of Gaussian eigenvalues
+//'
+//' Calculate eigenvalues of Hessian of Gaussian
+//' @param x Numeric matrix
+//' @param sigma Scale of smoothing
+//' @return List of 2 matrices with 1st and 2nd eigenvalue
+//' @export
+// [[Rcpp::export]]
+List  hessian_of_gaussian_eigenvalues(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());
+  MultiArray<2, TinyVector<float, 3> >  hessian(shape);
+  MultiArray<2, TinyVector<float, 2> >  eigen(shape);
+  hessianMatrixOfGaussian(input_image, hessian, sigma);
+  tensorEigenvaluesMultiArray(hessian, eigen);
+  MultiArrayView<2, float, StridedArrayTag> first = eigen.bindElementChannel(0);
+  MultiArrayView<2, float, StridedArrayTag> second = eigen.bindElementChannel(1);
+  NumericMatrix first_result(nrows, ncols);
+  NumericMatrix second_result(nrows, ncols);
+  std::copy(first.begin(), first.end(), first_result.begin());
+  std::copy(second.begin(), second.end(), second_result.begin());
+  return Rcpp::List::create(Rcpp::Named("ei1") = first_result,
+                                     Rcpp::Named("ei2") = second_result);
+}
diff --git a/src/laplacian_of_gaussian.cpp b/src/laplacian_of_gaussian.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f1e8600f9035754be22261d571ec211f324e5433
--- /dev/null
+++ b/src/laplacian_of_gaussian.cpp
@@ -0,0 +1,28 @@
+//[[Rcpp::depends(vigRa)]]
+#include <Rcpp.h>
+#include <vigra/multi_array.hxx>
+#include <vigra/convolution.hxx>
+using namespace Rcpp;
+using namespace vigra;
+
+
+//' Laplacian of Gaussian
+//'
+//' Calculates laplacian of Gaussian.
+//' @param x Numeric matrix
+//' @param sigma Scale of the Gaussian kernel
+//' @return Numeric matrix; smoothed version of x
+//' @export
+// [[Rcpp::export]]
+NumericMatrix laplacian_of_gaussian(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());
+  NumericMatrix result(nrows, ncols);
+  MultiArrayView<2, double> output_image(shape, result.begin());
+  laplacianOfGaussian(input_image, output_image, sigma);
+  return(result);
+}
+
+
diff --git a/src/structure_tensor_eigenvalues.cpp b/src/structure_tensor_eigenvalues.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..841d636823f30bc516ee498f28114787d680b574
--- /dev/null
+++ b/src/structure_tensor_eigenvalues.cpp
@@ -0,0 +1,50 @@
+//[[Rcpp::depends(vigRa)]]
+#include <Rcpp.h>
+#include <vigra/multi_array.hxx>
+#include <vigra/convolution.hxx>
+#include <vigra/tensorutilities.hxx>
+#include <vigra/mathutil.hxx>
+#include <vigra/multi_tensorutilities.hxx>
+using namespace Rcpp;
+using namespace vigra;
+
+//' Structure tensor eigenvalues (2 argument)
+//'
+//' Calculate eigenvalues of structure tensor
+//' @param x Numeric matrix
+//' @param inner_sigma Scale of inner smoothing
+//' @param outer_sigma Scale of outer smoothing
+//' @return List of 2 matrices with 1st and 2nd eigenvalue
+//' @export
+// [[Rcpp::export]]
+List  structure_tensor_eigenvalues_2(NumericMatrix x, double inner_sigma, double outer_sigma) {
+  const int nrows = x.nrow();
+  const int ncols = x.ncol();
+  Shape2 shape(nrows, ncols);
+  MultiArrayView<2, double> input_image(shape, x.begin());
+  MultiArray<2, TinyVector<float, 3> >  tensor(shape);
+  MultiArray<2, TinyVector<float, 2> >  eigen(shape);
+  structureTensor(input_image, tensor, inner_sigma, outer_sigma);
+  tensorEigenvaluesMultiArray(tensor, eigen);
+  MultiArrayView<2, float, StridedArrayTag> first = eigen.bindElementChannel(0);
+  MultiArrayView<2, float, StridedArrayTag> second = eigen.bindElementChannel(1);
+  NumericMatrix first_result(nrows, ncols);
+  NumericMatrix second_result(nrows, ncols);
+  std::copy(first.begin(), first.end(), first_result.begin());
+  std::copy(second.begin(), second.end(), second_result.begin());
+  return Rcpp::List::create(Rcpp::Named("ei1") = first_result,
+                                     Rcpp::Named("ei2") = second_result);
+}
+
+
+//' Structure tensor eigenvalues
+//'
+//' Calculate eigenvalues of structure tensor
+//' @param x Numeric matrix
+//' @param sigma Scale of inner smoothing; outer smoothing is sigma/2.
+//' @return List of 2 matrices with 1st and 2nd eigenvalue
+//' @export
+// [[Rcpp::export]]
+List  structure_tensor_eigenvalues(NumericMatrix x, double sigma) {
+  return structure_tensor_eigenvalues_2(x, sigma, sigma/2.0);
+}