# Cubature over Polygonal Domains

Numerical integration of continuously differentiable functions f(x,y) over simple closed polygonal domains. The following cubature methods are implemented: product Gauss cubature (Sommariva and Vianello, 2007, ), the simple two-dimensional midpoint rule (wrapping 'spatstat' functions), adaptive cubature for radially symmetric functions via line integrate() along the polygon boundary (Meyer and Held, 2014, , Supplement B), and integration of the bivariate Gaussian density based on polygon triangulation. For simple integration along the axes, the 'cubature' package is more appropriate.

The R package polyCub implements cubature (numerical integration) over polygonal domains. It solves the problem of integrating a continuously differentiable function f(x,y) over simple closed polygons.

For the special case where the domain is rectangular with sides parallel to the axes (such as a bounding box), the packages cubature and R2Cuba are more appropriate (cf. `CRAN Task View: Numerical Mathematics`).

## Installation

You can install polyCub from CRAN via:

To install the development version from the GitHub repository, use:

## Motivation

The polyCub package evolved from the need to evaluate integrals of so-called spatial interaction functions (e.g., a Gaussian or power-law kernel) over the observation region of a spatio-temporal point process (Meyer et al, 2012, Biometrics, https://doi.org/10.1111/j.1541-0420.2011.01684.x). Such an observation region is described by a polygonal boundary, representing, for example, the shape of a country or administrative district.

The integration task could be simplified by either assuming a trivial kernel, such as f(x,y)=1, or by simply replacing the polygonal with a rectangular domain, such as the bounding box of the polygon. However, these crude approximations can be avoided by using efficient numerical integration methods for polygonal domains:

1. Product Gauss cubature as proposed by Sommariva and Vianello (2007, BIT Numerical Mathematics, https://doi.org/10.1007/s10543-007-0131-2).

2. The simple two-dimensional midpoint rule via `as.im.function()` from the spatstat package.

3. For radially symmetric functions f(x,y) = f_r(||(x-x_0,y-y_0)||), numerical integration can be made much more efficient via line `integrate()` along the boundary of the polygonal domain (Meyer and Held, 2014, The Annals of Applied Statistics, https://doi.org/10.1214/14-AOAS743, Supplement B, Section 2.4).

4. For bivariate Gaussian densities, integrals over polygons can be solved accurately (but slowly) based on a triangulation of the domain (via `tristrip()` from the gpclib package) and combinations of Gaussian cumulative densities (via `pmvnorm()` from the mvtnorm package).

The dedicated R package polyCub was established in 2013 to provide implementations of the above cubature methods and facilitate their use in different projects. For example, polyCub powers epidemic models in surveillance and phylogeographic analyses in rase.

The four different cubature rules are exemplified below.

## Example

We consider a function f(x,y) which all of the above cubature methods can handle: an isotropic zero-mean Gaussian density. polyCub expects the function's implementation `f` to take a two-column coordinate matrix as its first argument (as opposed to separate arguments for the x and y coordinates):

We use a simple hexagon as polygonal integration domain, here specified via an `"xylist"` of vertex coordinates:

An image of the function and the integration domain can be produced using polyCub's rudimentary (but convenient) plotting utility:

#### Supported polygon representations

The integration domain is typically represented using a dedicated class for polygons, such as `"owin"` from package spatstat:

All of polyCub's cubature methods as well as `plotpolyf()` understand

• `"owin"` from spatstat,

• `"gpc.poly"` from package rgeos (or gpclib), and

• `"SpatialPolygons"` from package sp.

Note that the integration domain may consist of more than one polygon (including holes).

### 1. Product Gauss cubature: `polyCub.SV()`

The polyCub package provides an R-interfaced C-translation of "polygauss: Matlab code for Gauss-like cubature over polygons" (Sommariva and Vianello, 2013, http://www.math.unipd.it/~alvise/software.html). The cubature rule is based on Green's theorem and incorporates appropriately transformed weights and nodes of one-dimensional Gauss-Legendre quadrature in both dimensions, thus the name "product Gauss cubature". It is exact for all bivariate polynomials if the number of cubature nodes is sufficiently large (depending on the degree of the polynomial).

For the above example, a reasonable approximation is already obtained with degree `nGQ = 3` of the one-dimensional Gauss-Legendre quadrature:

The involved nodes (displayed in the figure above) and weights can be extracted by calling `polyCub.SV()` with `f = NULL`, e.g., to determine the number of nodes:

For illustration, we create a variant of `polyCub.SV()`, which returns the number of function evaluations as an attribute:

We can use this function to investigate how the accuracy of the approximation depends on the degree `nGQ` and the associated number of cubature nodes:

### 2. Two-dimensional midpoint rule: `polyCub.midpoint()`

The two-dimensional midpoint rule in polyCub is a simple wrapper around `as.im.function()` and `integral.im()` from package spatstat. The polygon is converted to a binary pixel image and the integral is approximated as the sum of (pixel area * f(pixel midpoint)) over all pixels whose midpoint is part of the polygon.

Using a pixel size of `eps = 0.5` (here yielding 270 pixels), we obtain:

### 3. Adaptive cubature for isotropic functions: `polyCub.iso()`

A radially symmetric function can be expressed in terms of the distance r from its point of symmetry: f(r). If the antiderivative of r times f(r), called `intrfr()`, is analytically available, Green's theorem leads us to a cubature rule which only needs one-dimensional numerical integration. More specifically, `intrfr()` will be `integrate()`d along the edges of the polygon. The mathematical details are given in https://doi.org/10.1214/14-AOAS743SUPPB (Section 2.4).

For the bivariate Gaussian density `f` defined above, the integral from 0 to R of `r*f(r)` is analytically available as:

With this information, we can apply the cubature rule as follows:

Note that we do not even need the original function `f`.

If `intrfr()` is missing, it can be approximated numerically using `integrate()` for `r*f(r)` as well, but the overall integration will then be much less efficient than product Gauss cubature.

Package polyCub exposes a C-version of `polyCub.iso()` for use by other R packages (notably surveillance) via `LinkingTo: polyCub` and `#include <polyCubAPI.h>`. This requires the `intrfr()` function to be implemented in C as well. See https://github.com/bastistician/polyCub/blob/master/tests/testthat/polyiso_powerlaw.c for an example.

### 4. Integration of the bivariate Gaussian density: `polyCub.exact.Gauss()`

Abramowitz and Stegun (1972, Section 26.9, Example 9) offer a formula for the integral of the bivariate Gaussian density over a triangle with one vertex at the origin. This formula can be used after triangulation of the polygonal domain via `tripstrip()` from the gpclib package. The core of the formula is an integral of the bivariate Gaussian density with zero mean, unit variance and some correlation over an infinite rectangle [h, Inf] x [0, Inf], which can be computed accurately using `pmvnorm()` from the mvtnorm package.

For the above example, we obtain:

The required triangulation as well as the numerous calls of `pmvnorm()` make this integration algorithm quiet cumbersome. For large-scale integration tasks, it is thus advisable to resort to the general-purpose product Gauss cubature rule `polyCub.SV()`.

Note: There is also a function `circleCub.Gauss()` to calculate the integral of an isotropic Gaussian density over a circular domain (which requires nothing more than a single call of `pchisq()`).

## Feedback

Contributions are welcome! Please submit suggestions or report bugs at https://github.com/bastistician/polyCub/issues. You can also reach me via e-mail ([email protected]). Note that pull requests should only be submitted after a discussion of the underlying issue.

The polyCub package is free and open source software, licensed under the GPLv2.

# polyCub 0.7.0 (2018-10-11)

• Package polyCub no longer attaches package sp (moved from "Depends" to "Imports").

• The R code of the examples is no longer installed by default. Use the `--example` flag of R CMD INSTALL to achieve that.

• The README now exemplifies the four different cubature rules.

# polyCub 0.6.1 (2017-10-02)

• The exported C-function `polyCub_iso()` ...

• did not handle its `stop_on_error` argument correctly (it would always stop on error).

• now detects non-finite `intrfr` function values and gives an informative error message (rather than just reporting "abnormal termination of integration routine").

• Package polyCub no longer strictly depends on package spatstat. It is only required for `polyCub.midpoint()` and for polygon input of class `"owin"`.

# polyCub 0.6.0 (2017-05-24)

• Added full C-implementation of `polyCub.iso()`, which is exposed as `"polyCub_iso"` for use by other R packages (notably future versions of surveillance) via `LinkingTo: polyCub` and `#include <polyCubAPI.h>`.

• Accommodate CRAN checks: add missing import from graphics, register native routines and disable symbol search

# polyCub 0.5-2 (2015-02-25)

• `polyCub.midpoint()` works directly with input polygons of classes `"gpc.poly"` and `"SpatialPolygons"`, since package polyCub now registers corresponding `as.owin`-methods.

• `polyCub.exact.Gauss()` did not work if the `tristrip` of the transformed input polygon contained degenerate triangles (spotted by Ignacio Quintero).

• Line integration in `polyCub.iso()` could break due to division by zero if the `center` point was part of the polygon boundary.

# polyCub 0.5-1 (2014-10-24)

• Nodes and weights for `polyCub.SV()` were only cached up to `nGQ=59`, not 60 as announced in version 0.5-0. Fixed that which also makes examples truly run without statmod.

• In `polyCub.SV()`, the new special setting `f=NULL` means to only compute nodes and weights.

• Internal changes to the `"gpc.poly"` converters to accommodate spatstat 1.39-0.

# polyCub 0.5-0 (2014-05-07)

• `polyCub.SV()` gained an argument `engine` to choose among available implementations. The new and faster C-implementation is the default. There should not be any numerical differences in the result of the cubature.

• Package statmod is no longer strictly required (imported). Nodes and weights for Gauss-Legendre quadrature in `polyCub.SV()` are now cached in the polyCub package up to `nGQ=60`. statmod`::gauss.quad` is only queried for a higher number of nodes.

# polyCub 0.4-3 (2014-03-14)

• `polyCub.iso()` ...

• could not handle additional arguments for `integrate()` given in the `control` list.

• now applies the `control` arguments also to the numerical approximation of `intrfr`.

• The `checkintrfr()` function is exported and documented.

# polyCub 0.4-2 (2014-02-12)

• `plotpolyf()` ...

• gained an additional argument `print.args`, an optional list of arguments passed to `print.trellis()` if `use.lattice=TRUE`.

• passed a data frame of coordinates to `f` instead of a matrix as documented.

# polyCub 0.4-1 (2013-12-05)

• This version solely fixes a missing NAMESPACE import to make package polyCub again compatible with older versions of spatstat (< 1.33-0).

# polyCub 0.4-0 (2013-11-19)

## INFRASTRUCTURE

• rgeos (and therefore the GEOS library) is no longer strictly required (moved from "Imports" to "Suggests").

• Added `coerce`-methods from `"Polygons"` (or `"SpatialPolygons"` or `"Polygon"`) to `"owin"` (`as(..., "owin")`).

• S4-style `coerce`-methods between `"gpc.poly"` and `"Polygons"`/`"owin"` have been removed from the package (since we no longer import the formal class `"gpc.poly"` from gpclib or rgeos). However, there are two new functions `gpc2owin` and `owin2gpc` similar to those dropped from spatstat since version 1.34-0.

• Moved `discpoly()` back to surveillance since it is only used there.

• The latter two changes cause surveillance version 1.6-0 to be incompatible with this new version of polyCub. Appropriate modifications have been made in the new version 1.7-0 of surveillance.

## SPEED-UP `polyCub.SV()`

• thorough optimization of `polyCub.SV()`-related code resulted in about 27% speed-up:

• use `mapply()` instead of a `for`-loop

• avoid `cbind()`

• use `tcrossprod()`

• less object copying

## MINOR CHANGES

• `xylist()` is now exported. It simply extracts polygon coordinates from various spatial classes (with same unifying intention as `xy.coords()`).

• A `polyregion` of class `"SpatialPolygons"` of length more than 1 now works in `polyCub`-methods.

• Use aspect ratio of 1 in `plotpolyf()`.

# polyCub 0.3-1 (2013-08-22)

• This version solely fixes a few typos and a technical note from `R CMD check` in the current R development version (also import packages into the NAMESPACE which are listed in the "Depends" field).

# polyCub 0.3-0 (2013-07-06)

• New cubature method `polyCub.iso()` specific to isotropic functions (thanks to Emil Hedevang for the basic idea).

• New function `plotpolyf()` to plot a polygonal domain on top of an image of a bivariate function.

• The package now depends on R >= 2.15.0 (for `.rowSums()`).

• The package no longer registers `"owin"` as an S4-class since we depend on the sp package which does the job. This avoids a spurious warning (in `.simpleDuplicateClass()`) upon package installation.

• In `discpoly()`, the argument `r` has been renamed to `radius`. This is backward compatible by partial argument matching in old code.

# polyCub 0.2-0 (2013-05-09)

• This is the initial version of the polyCub package mainly built on functions previously maintained within the surveillance package. These methods for cubature of polygonal domains have been outsourced into this separate polyCub package since they are of general use for other packages as well.

• The polyCub package has more documentation and tests, avoids the use of gpclib as far as possible (using rgeos instead), and solves a compatibility issue with package maptools (use `setClass("owin")` instead of `setOldClass("owin")`).

# Reference manual

install.packages("polyCub")

0.7.1 by Sebastian Meyer, 11 days ago

https://github.com/bastistician/polyCub

Report a bug at https://github.com/bastistician/polyCub/issues

Browse source code at https://github.com/cran/polyCub

Authors: Sebastian Meyer [aut, cre, trl] , Leonhard Held [ths] , Michael Hoehle [ths]

Documentation:   PDF Manual