Rcpp – Integrating C++ Into R

Using Rcpp

Following my last post I will now show you how to integrate C++ code into a R package using Rcpp. Nothing revolutionary to do, but maybe it will still help some of you. 🙂
I will demonstrate this on my GitHub repository and R package BiocStyle::Githubpkg("David-J-R/Raspository"). And I will progress on optimizing the Dithering method I implemented last time. For the following I suppose, that you already use roxygen2.

Motivation

So at first you should ask yourself why you wanna introduce C++ code to your package. One very valid reason might be, that the C++ code is already there and you just want to integrate it into R or maybe even just write a wrapper in R around it.
The main reason why I did this was to speed up my code. Some operations in R are already pretty fast. But if you want to write a function with a loop that does complicated things, that would be torture to implement with apply functions, if even possible, C++ is a good choice.

There might of course be even more reasons, but we will save them for another time…
Let’s now get to how to do it.

Preparation For Integrating Rcpp Into R Packages

Before we can seamlessly integrate C++ using into our package we need to do some preparations. First we have to import the Rcpp package into our namespace and link to its dynamic library.
For this cause we add the following to the DESCRIPTION file of the package:

Imports: Rdpack, BEclear, futile.logger, stats, dequer, jpeg,    data.table, RobustRankAggreg, igraph, graphics, grDevices, stats, methods, abind, Rcpp
LinkingTo: Rcpp

And we also have to use the dynamic that will be created in our package. For this cause we add the following statement to the Raspository.R:

#' Raspository-package
#' 
#' @author David Rasp
#' 
#' @docType package
#' 
#' @import Rdpack
#' 
#' @title The Raspository
#' 
#' @useDynLib Raspository
"_PACKAGE"

As far as I understand it you could also add this statement anywhere into the documentation of your package.

Creating our function with Rcpp

So the translation of the dissipatePixel function from last time is pretty straight forward. It’s just changing the syntax a bit. I also added a roxygen skeleton. It also works the same as in R, just that you use two backslashes instead of a hash-tag to start a line. The namespace of Rcpp provides you with classes like NumericMatrix, that can be used similar to their R counterparts and the objects in R can be treated as them. If you want some more references, I would advise you to read the Rcpp quickref by Dirk Eddelbuettel. I have it open most of the time I’m working with it.

#include <Rcpp.h>
using namespace Rcpp;


//' Calculating the Dissipation of Pixels
//' 
//' @export dissipatePixel
//' @import Rcpp
//' 
//' 
//' @references \insertRef{Hagenburg2009}{Raspository}
//' 
//' @return the image as matrix after the next time step
//'
// [[Rcpp::export]]
SEXP dissipatePixel(const SEXP& imgOriginal, const double& minimalTreshold) {
    NumericMatrix img(imgOriginal);
    NumericMatrix imgAtNewStep(img.nrow(), img.ncol());

    for(std::size_t i = 0; i < img.nrow(); i++){
        for(std::size_t j = 0; j < img.ncol(); j++){

            double currentPixel = img(i,j);
            double dissipated = 0.0;

            if(currentPixel > 1.0 || currentPixel < minimalTreshold){

                double toDissipate = currentPixel; 

                if(currentPixel > 1){
                    toDissipate -= 1.0;
                }

                double tmp1 = toDissipate / 9.0;
                double tmp2 = toDissipate / 36.0;

                if(i > 0){
                    imgAtNewStep(i - 1, j) += tmp1;
                    dissipated += tmp1;
                }

                if(j > 0){
                    imgAtNewStep(i,j - 1) += tmp1;
                    dissipated += tmp1;
                }

                if(i < img.nrow() - 1){
                    imgAtNewStep(i + 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if(j < img.ncol() - 1){
                    imgAtNewStep(i,j + 1) += tmp1;
                    dissipated += tmp1;
                }

                if( i > 0 && j > 0){
                    imgAtNewStep(i - 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i > 0 && j < img.ncol() - 1){
                    imgAtNewStep(i - 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > 0){
                    imgAtNewStep(i + 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > img.ncol() - 1){
                    imgAtNewStep(i + 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }


            }else{

                double tmp1 = currentPixel / 9.0;
                double tmp2 = currentPixel / 36.0;

                if( i > 1 && img(i - 1,j) > currentPixel && img(i - 1,j) < 1){
                    imgAtNewStep(i - 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if( j > 0 && img(i,j - 1) > currentPixel && img(i,j - 1) < 1){
                    imgAtNewStep(i,j - 1) += tmp1;
                    dissipated += tmp1;
                }

                if(i < img.nrow() - 1 && img(i + 1,j) > currentPixel && img(i + 1,j) < 1){
                    imgAtNewStep(i + 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if(j < img.ncol() - 1 && img(i,j + 1) > currentPixel && img(i,j + 1) < 1){
                    imgAtNewStep(i,j + 1) += tmp1;
                    dissipated += tmp1;
                }

                if( i > 0 && j > 0 && img(i - 1,j - 1) > currentPixel && img(i - 1,j - 1) < 1){
                    imgAtNewStep(i - 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i > 0 && j < img.ncol() - 1 && img(i - 1,j + 1) > currentPixel && img(i - 1,j + 1) < 1){
                    imgAtNewStep(i - 1,j + 1) += 1.0 / 36.0 * currentPixel;
                    dissipated += 1.0 / 36.0 * currentPixel;
                }

                if( i < img.nrow() - 1 && j > 0 && img(i + 1,j - 1) > currentPixel && img(i + 1,j - 1) < 1){
                    imgAtNewStep(i + 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > img.ncol() - 1 && img(i + 1,j + 1) > currentPixel && img(i + 1,j + 1) < 1){
                    imgAtNewStep(i + 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }

            }

            // add the non dissipated amount to the same pixel in next time-step
            imgAtNewStep(i,j) += currentPixel - dissipated;
        }
    }

return imgAtNewStep;
}

After writing and saving this function, we have to generate the RccpExports.cpp RccpExports.R files, which basically contains a wrapper for our C++ function, which then after compiling the code can be loaded in R.
We can comfortably do this with:

Rcpp::compileAttributes()

Then we should roxygenation our package, so that the namespace and the documentation get updated. We do this with:

roxygen2::roxygenise()

Now finally we can recompile the whole package and load the new version. The simplest way imo is to just use Ctrl+Shift+B in RStudio.

Performance comparison

Of we should also check, how much faster the new version is. I’ll use microbenchmark for this cause.

library(microbenchmark)
data <- matrix(runif(625), 25, 25)
microbenchmark(dissipatePixel(data, minimalTreshold = 0.05), dissipatePixelOld(data, minimalTreshold = 0.05))
## Unit: microseconds
##                                             expr      min       lq
##     dissipatePixel(data, minimalTreshold = 0.05)   70.891   72.368
##  dissipatePixelOld(data, minimalTreshold = 0.05) 4479.941 4622.600
##        mean   median        uq        max neval
##    95.23017   77.003   82.4445   1486.321   100
##  6514.79719 4677.444 4779.1125 161310.562   100

I think the result speaks for itself. This optimization was very worth while. The version implemented in C++ is about 60 times faster. That’s a lot. This could be the difference between one second and a minute. And of course longer calculations also suck more electricity, which means they have a larger CO2 fingerprint. So be nice to our atmosphere and code effective. 😉

That’s it for now.
Later,
David.

Availability Of The Code

You can access a maintained version of the code for the color spaces in my GitHub repository Raspository under src/dissipatePixel.cpp.

Please follow and like us:
error0

Leave a Reply

Your email address will not be published. Required fields are marked *