ppdist: Compute Distance Between Two Point Patterns

Description

Based on an arbitrary matrix of "distances" between the points of two point patterns $\xi$ and $\eta$, this function computes versions of the transport-transform distance between $\xi$ and $\eta$.

Usage

ppdist(
  dmat,
  penalty = 1,
  type = c("tt", "rtt", "TT", "RTT"),
  ret_matching = FALSE,
  p = 1,
  precision = NULL
)

Value

The corresponding distance between the point patterns if ret_matching

is FALSE.

Otherwise a list with components dist containing this distance and two vectors target1, target2 of integers, where

target$i$ specifies the indices of the points in the other pattern that the points of the $i$-th pattern are matched to and NA every time a point is deleted.

There may be a minus in front of an index, where

-j indicates that the corresponding pairing with point j

would be over a distance of more than $2^{1/p} \cdot \code{penalty}$. This is equivalent to saying that the corresponding point of the first pattern is deleted and the $j$-th point of the second pattern is added.

Note that having more than one minus implies that the matching is non-unique.

Arguments

dmat: a matrix specifying in its $(i,j)$-th entry the distance from the i-th point of $\xi$ to the j-th point of $\eta$.
penalty: a positive number. The penalty for adding/deleting points.
type: either "tt"/"TT" for the transport-transform metric or "rtt"/"RTT" for the relative transport-transform metric.
ret_matching: logical. Shall the optimal point matching be returned?
p: a number $>0$. The matching is chosen such that the p-th order sum ($\ell_p$-norm) is minimized.
precision: a small positive integer value. The precisions of the computations, which are currently performed in integers. After correcting for the penalty, dmat^p is divided by its largest entry, multiplied by 10^precision and rounded to compute the optimal matching. The default value NULL chooses maximal integer precision possible, which is precision = 9 on almost all systems.

Author

Dominic Schuhmacher schuhmacher@math.uni-goettingen.de

Details

The transport-transform (TT) distance gives the minimal total cost for “morphing” the pattern $\xi$ into the pattern $\eta$ by way of shifting points (at costs specified in dmat) and adding or deleting points (each at cost penalty). The total cost is determined as $$\biggl( \sum_{j=1}^n c_j^p \biggr)^{1/p},$$ where $c_j$ denotes the cost for the $j$th individual operation and $n$ is the cardinality of the larger point pattern.

The relative transport-transform (RTT) metric is exactly the same, but the sum in the total cost is divided by the larger cardinality: $$\biggl( \frac{1}{n}\sum_{j=1}^n c_j^p \biggr)^{1/p}.$$

The TT- and RTT-metrics form an umbrella concept that includes the OSPA and Spike Time metrics frequently used in the literature. See Müller, Schuhmacher and Mateu (2020) for details.

References

Raoul Müller, Dominic Schuhmacher and Jorge Mateu (2020).
Metrics and Barycenters for Point Pattern Data.
Statistics and Computing 30, 953-972.
tools:::Rd_expr_doi("10.1007/s11222-020-09932-y")

Examples

Run this code

  # small example
  # -------------
  set.seed(181230)
  xi <- spatstat.random::rpoispp(20)
  eta <- spatstat.random::rpoispp(20)
  dmat <- spatstat.geom::crossdist(xi,eta)
  res <- ppdist(dmat, penalty=1,  type="rtt", ret_matching=TRUE, p=1)
  plotmatch(xi, eta, dmat, res, penalty=1, p=1)
  res$dist

  # for comparison: ospa-distance computation from spatstat:
  res_ospa <- spatstat.geom::pppdist(xi,eta,"spa")
  res_ospa$distance  # exactly the same as above because nothing gets cut off 


  # same example, but with a much smaller penalty for adding/deleting points
  # --------------------------------------------------------------- 
  res <- ppdist(dmat, penalty=0.1,  type="rtt", ret_matching=TRUE, p=1)
  plotmatch(xi, eta, dmat, res, penalty=0.1, p=1)
    # dashed lines indicate points that are deleted and re-added at new position
    # grey segments on dashed lines indicate the cost of deletion plus re-addition
  res$dist
  
  # for comparison: ospa-distance computation from spatstat
  # (if things do get cut off, we have to ensure that the cutoff distances
  # are the same, thus cutoff = 2^(1/p) * penalty):
  res_ospa <- spatstat.geom::pppdist(xi,eta,"spa",cutoff=0.2)
  res_ospa$distance  # NOT the same as above
  res_ospa$distance - abs(xi$n-eta$n) * 0.1 / max(xi$n,eta$n)  # the same as above
  
  
  # a larger example
  # --------------------------------------------------------------- 
  set.seed(190203)
  xi <- spatstat.random::rpoispp(2000)
  eta <- spatstat.random::rpoispp(2000)
  dmat <- spatstat.geom::crossdist(xi,eta)
  res <- ppdist(dmat, penalty = 0.1,  type = "rtt", ret_matching = TRUE, p = 1)
  res$dist
  # takes about 2-3 seconds

Run the code above in your browser using DataLab