polyCub (version 0.7.1)

polyCub.midpoint: Two-Dimensional Midpoint Rule

Description

The surface is converted to a binary pixel image using the as.im.function method from package spatstat (Baddeley and Turner, 2005). The integral under the surface is then approximated as the sum over (pixel area * f(pixel midpoint)).

Usage

polyCub.midpoint(polyregion, f, ..., eps = NULL, dimyx = NULL,
  plot = FALSE)

Arguments

polyregion

a polygonal integration domain. It can be any object coercible to the spatstat class "owin" via a corresponding as.owin-method. Note that this includes polygons of the classes "gpc.poly" and "'>SpatialPolygons", because polyCub defines methods as.owin.gpc.poly and as.owin.SpatialPolygons, respectively.

f

a two-dimensional real-valued function. As its first argument it must take a coordinate matrix, i.e., a numeric matrix with two columns, and it must return a numeric vector of length the number of coordinates.

...

further arguments for f.

eps

width and height of the pixels (squares), see as.mask.

dimyx

number of subdivisions in each dimension, see as.mask.

plot

logical indicating if an illustrative plot of the numerical integration should be produced.

Value

The approximated value of the integral of f over polyregion.

References

Baddeley, A. and Turner, R. (2005). spatstat: an R package for analyzing spatial point patterns. Journal of Statistical Software, 12 (6), 1-42.

See Also

Other polyCub-methods: polyCub.SV, polyCub.exact.Gauss, polyCub.iso, polyCub

Examples

# NOT RUN {
## a function to integrate (here: isotropic zero-mean Gaussian density)
f <- function (s, sigma = 5)
    exp(-rowSums(s^2)/2/sigma^2) / (2*pi*sigma^2)

## a simple polygon as integration domain
hexagon <- list(
    list(x = c(7.33, 7.33, 3, -1.33, -1.33, 3),
         y = c(-0.5, 4.5, 7, 4.5, -0.5, -3))
)

if (require("spatstat")) {
    hexagon.owin <- owin(poly = hexagon)

    show_midpoint <- function (eps)
    {
        plotpolyf(hexagon.owin, f, xlim = c(-8,8), ylim = c(-8,8),
                  use.lattice = FALSE)
        ## add evaluation points to plot
        with(as.mask(hexagon.owin, eps = eps),
             points(expand.grid(xcol, yrow), col = t(m), pch = 20))
        title(main = paste("2D midpoint rule with eps =", eps))
    }

    ## show nodes (eps = 0.5)
    show_midpoint(0.5)

    ## show pixel image (eps = 0.5)
    polyCub.midpoint(hexagon.owin, f, eps = 0.5, plot = TRUE)

    ## use a decreasing pixel size (increasing number of nodes)
    for (eps in c(5, 3, 1, 0.5, 0.3, 0.1))
        cat(sprintf("eps = %.1f: %.7f\n", eps,
                    polyCub.midpoint(hexagon.owin, f, eps = eps)))
}
# }