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.geom (Baddeley et al., 2015). 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)

Value

The approximated value of the integral of f over polyregion.

Arguments

polyregion: a polygonal integration domain. It can be any object coercible to the spatstat.geom 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. sf also registers suitable as.owin methods for its "(MULTI)POLYGON" classes.
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.

References

Baddeley A, Rubak E, Turner R (2015). Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press, London.

Examples

Run this code

## 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.geom")) {
    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)))
}