ppdist {ttbary} | R Documentation |
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
)
Arguments
dmat |
a matrix specifying in its |
penalty |
a positive number. The penalty for adding/deleting points. |
type |
either |
ret_matching |
logical. Shall the optimal point matching be returned? |
p |
a number |
precision |
a small positive integer value. The precisions of the computations, which
are currently performed in integers. After correcting for the penalty, |
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.
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.
Author(s)
Dominic Schuhmacher schuhmacher@math.uni-goettingen.de
References
Raoul Müller, Dominic Schuhmacher and Jorge Mateu (2020).
Metrics and Barycenters for Point Pattern Data.
Statistics and Computing 30, 953-972.
doi:10.1007/s11222-020-09932-y
Examples
# 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