Learn R Programming

polyCub (version 0.9.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.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.

See Also

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

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)))
}

Run the code above in your browser using DataLab