Add support for threads argument and related option

DESCRIPTION

@@ -16,4 +16,4 @@ License: GPL-2
 URL: https://github.com/jonclayden/mmand
 BugReports: https://github.com/jonclayden/mmand/issues
 Roxygen: list(old_usage=TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3

R/distance.R

@@ -19,6 +19,8 @@
 #'   is returned, such that distances from the region boundary are negative
 #'   within the region and positive outside. Otherwise, distances are zero
 #'   within the region.
+#' @param threads If a positive integer, and the package is compiled with
+#'   OpenMP support, the number of threads to use during the calculation.
 #' @return An array of the same dimension as the original, whose elements give
 #'   the Euclidean distance from that element to the nearest "on" element in
 #'   the original.
@@ -41,7 +43,7 @@ distanceTransform <- function (x, ...)
 
 #' @rdname distanceTransform
 #' @export
-distanceTransform.default <- function (x, pixdim = TRUE, signed = FALSE, ...)
+distanceTransform.default <- function (x, pixdim = TRUE, signed = FALSE, threads = getOption("mmand.threads"), ...)
 {
     x <- as.array(x)
     if (!is.numeric(x) && !is.logical(x))
@@ -71,10 +73,10 @@ distanceTransform.default <- function (x, pixdim = TRUE, signed = FALSE, ...)
         }
         else
             value <- attr(isBinary, "value")
-        returnValue <- .Call(C_distance_transform, x, pixdim) - .Call(C_distance_transform, value-x, pixdim)
+        returnValue <- .Call(C_distance_transform, x, pixdim, threads) - .Call(C_distance_transform, value-x, pixdim, threads)
     }
     else
-        returnValue <- .Call(C_distance_transform, x, pixdim)
+        returnValue <- .Call(C_distance_transform, x, pixdim, threads)
 
     if (length(dim(x)) > 1)
         dim(returnValue) <- dim(x)

R/resample.R

@@ -14,6 +14,8 @@
 #'   each resampled value, or the name of one.
 #' @param pointType A string giving the type of the point specification being
 #'   used. Usually can be left as \code{"auto"}.
+#' @param threads If a positive integer, and the package is compiled with
+#'   OpenMP support, the number of threads to use during the calculation.
 #' @param factor A vector of scale factors, which will be recycled to the
 #'   dimensionality of \code{x}.
 #' @param \dots Additional options, such as kernel parameters.
@@ -34,7 +36,7 @@ resample <- function (x, points, kernel, ...)
 
 #' @rdname resample
 #' @export
-resample.default <- function (x, points, kernel, pointType = c("auto","general","grid"), ...)
+resample.default <- function (x, points, kernel, pointType = c("auto","general","grid"), threads = getOption("mmand.threads"), ...)
 {
     x <- as.array(x)
     if (!is.numeric(x) && !is.logical(x))
@@ -68,7 +70,7 @@ resample.default <- function (x, points, kernel, pointType = c("auto","general",
     else if (is.list(points))
         points <- lapply(points, "-", 1)
 
-    result <- .Call(C_resample, x, kernel, list(type=pointType,points=points))
+    result <- .Call(C_resample, x, kernel, list(type=pointType,points=points), threads)
 
     if (is.list(points) && nDims > 1)
         dim(result) <- sapply(points, length)

src/main.cpp

@@ -5,6 +5,10 @@
 #include "Resampler.h"
 #include "Morpher.h"
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 using namespace Rcpp;
 using namespace std;
 
@@ -127,7 +131,7 @@ BEGIN_RCPP
 END_RCPP
 }
 
-RcppExport SEXP resample (SEXP data_, SEXP kernel_, SEXP samplingScheme_)
+RcppExport SEXP resample (SEXP data_, SEXP kernel_, SEXP samplingScheme_, SEXP threads_)
 {
 BEGIN_RCPP
     Array<double> *array = arrayFromData(data_);
@@ -137,6 +141,11 @@ BEGIN_RCPP
     List samplingScheme(samplingScheme_);
     string schemeType = as<string>(samplingScheme["type"]);
 
+#ifdef _OPENMP
+    if (!Rf_isNull(threads_) && as<int>(threads_) > 0)
+        omp_set_num_threads(as<int>(threads_));
+#endif
+
     if (schemeType.compare("general") == 0)
     {
         NumericMatrix points = samplingScheme["points"];
@@ -228,11 +237,15 @@ BEGIN_RCPP
 END_RCPP
 }
 
-RcppExport SEXP distance_transform (SEXP data_, SEXP _usePixdim)
+RcppExport SEXP distance_transform (SEXP data_, SEXP usePixdim_, SEXP threads_)
 {
 BEGIN_RCPP
     Array<double> *array = arrayFromData(data_);
-    Distancer distancer(array, as<bool>(_usePixdim));
+#ifdef _OPENMP
+    if (!Rf_isNull(threads_) && as<int>(threads_) > 0)
+        omp_set_num_threads(as<int>(threads_));
+#endif
+    Distancer distancer(array, as<bool>(usePixdim_));
     Array<double> *distances = distancer.run();
     SEXP result = wrap(distances->getData());
     delete distances;
@@ -245,10 +258,10 @@ static R_CallMethodDef callMethods[] = {
     { "is_symmetric",           (DL_FUNC) &is_symmetric,            1 },
     { "get_neighbourhood",      (DL_FUNC) &get_neighbourhood,       2 },
     { "sample_kernel",          (DL_FUNC) &sample_kernel,           2 },
-    { "resample",               (DL_FUNC) &resample,                3 },
+    { "resample",               (DL_FUNC) &resample,                4 },
     { "morph",                  (DL_FUNC) &morph,                   6 },
     { "connected_components",   (DL_FUNC) &connected_components,    2 },
-    { "distance_transform",     (DL_FUNC) &distance_transform,      2 },
+    { "distance_transform",     (DL_FUNC) &distance_transform,      3 },
     { NULL, NULL, 0 }
 };