[Rd] inline C/C++ in R: question and suggestion

From: Oleg Sklyar <osklyar_at_ebi.ac.uk>
Date: Tue, 22 May 2007 18:59:00 +0100


This is a question and maybe an announcement.

We've been discussing in the group that it would be nice to have a mechanism for something like "inline" C/C++ function calls in R. I do not want to reinvent the wheel, therefore, if something like that already exists, please give me a hint -- I could not find anything. If not, here is a working solution, please criticise so I could improve it.

Example: I work on images (Bioconductor:EBImage) and just came to a point where I need to apply certain functions to image data, which are grey scale intensities in the range [0,1]. Assume I want to transform my image data from I(x,y,i) to exp(-(d/s)^2)*I(x,y,i), where I is the original intensity in dependence on coordinates x,y and frame i; s is a given value and d^2=(x-centre.x)^2+(y-centre.y)^2 for a given centre. Trying an R loop will run forever already on moderate image sizes as I do not see how to vectorize it.

Now, below is the solution using the "inline" C code, completely in R and runs instantly. I created a small package "inline" that simply encapsulates a quite simple function "cfunction". The package source is available from http://www.ebi.ac.uk/~osklyar/inline -- please give it a try and I would be happy to hear your comments, both on already existing implementations for "inline" calls and on the current one. It is a draft and I will be happy to improve it (especially comments on how to ensure that the shared object is unloaded when the function is removed). In the example below EBImage is not required, I use randomly generated values instead of images, but the output it quite obvious. After installing "inline" the example should just work by copy-pasting.

Best and thanks in advance,
Oleg

code <- character(17)

code[1] <- "  SEXP res;"
code[2] <- "  int nprotect = 0, nx, ny, nz, x, y;"
code[3] <- "  PROTECT(res = Rf_duplicate(a)); nprotect++;"
code[4] <- "  nx = INTEGER(GET_DIM(a))[0];"
code[5] <- "  ny = INTEGER(GET_DIM(a))[1];"
code[6] <- "  nz = INTEGER(GET_DIM(a))[2];"
code[7] <- "  double sigma2 = REAL(s)[0] * REAL(s)[0], d2 ;"
code[8] <- "  double cx = REAL(centre)[0], cy = REAL(centre)[1], *data,
*rdata;"
code[9] <- " for (int im = 0; im < nz; im++) {" code[10] <- " data = &(REAL(a)[im*nx*ny]); rdata = &(REAL(res)[im*nx*ny]);"
code[11] <- "    for (x = 0; x < nx; x++)"
code[12] <- "       for (y = 0; y < ny; y++) {"
code[13] <- "         d2 = (x-cx)*(x-cx) + (y-cy)*(y-cy);"
code[14] <- "         rdata[x + y*nx] = data[x + y*nx] * exp(-d2/sigma2);"
code[15] <- "  }}"
code[16] <- "  UNPROTECT(nprotect);"
code[17] <- "  return res;"

library(inline)

funx <- cfunction(signature(a="array", s="numeric", centre="numeric"), code)

x <- array(runif(50*50), c(50,50,1))

res <- funx(a=x, s=15, centre=c(30,35))
image(res[,,1])
res <- funx(x, 10, c(15,15))
x11(); image(res[,,1])

-- 
Dr Oleg Sklyar * EBI/EMBL, Cambridge CB10 1SD, England * +44-1223-494466

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Tue 22 May 2007 - 18:00:50 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 22 May 2007 - 20:35:46 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.