On 23/03/15 23:44, lluis.hurt...@uv.es wrote:
Dear all, I am currently working with the spatstat package with 3D samples. I am trying to evaluate a non analytical function over the window that encloses the sample and I need to know which is the fastest way of doing it.
I gather that you mean that you want to evaluate the *integral* of the function over the window.
What do you mean by "non analytical"? You have, it would appear, some way of calculating the function at an arbitrary point (x,y,z) in 3-space, which makes it at least in some some sense "analytical".
The function input is a 3 coordinate position in the window (x,y,z) and a list of parameters (a,b,c). The output is a numerical value. n <- function(x,y,z,a,b,c) But I need to do it over the whole volume. For 2 dimensions it can be done with A <- as.im(function,window,parameters) norm <- integral.im(A).
There is nothing very sophisticated about the way integral.im() works. Just look at the code. The operative line of code is:
a <- with(f, sum(v, na.rm = TRUE) * xstep * ystep)
For 3 dimensions I have tried to pass an array of a grid covering the window (like a quadrature scheme) and then summing up the output array, but I would like to know if there is any faster way of integrating the function.
It sounds, at least vaguely, as though you are imitating in 3D what integral.im() does in 2D. So, given that you are doing it right, you are doing as well as integral.im() does.
The calculations should be pretty fast since sum() is pretty fast. The bulk of the computational burden is likely to be the evaluation of the function at all of the voxel centres.
Note that this is a very unsophisticated method of numerical quadrature. The integrand is being approximated by a step function. The method seems to work well enough for the uses that we put it to in spatstat. Whether it works well enough for your purposes depends on what your purposes are.
cheers, Rolf Turner -- Rolf Turner Technical Editor ANZJS Department of Statistics University of Auckland Phone: +64-9-373-7599 ext. 88276 Home phone: +64-9-480-4619 ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
