How to measure similarity of SpatialLines objects

I created two `SpatialLines` objects in R: .

These objects were created this way:

``````library(sp)
xy <- cbind(x,y)
xy.sp = sp::SpatialPoints(xy)
spl1 <- sp::SpatialLines(list(Lines(Line(xy.sp), ID="a")))
``````

Now I want to somehow conclude that this is the same line rotated and flipped, and that the difference between them is equal to 0 (i.e. shape is equal).

To do that, one can use `maptools` package and rotate line #1, e.g.:

``````spl180 <- maptools::elide(spl1, rotate=180)
``````

Each rotated line must then be checked versus line #2 using `rgeos` package, e.g.:

``````hdist <- rgeos::gDistance(spl180, spl2, byid=FALSE, hausdorff=TRUE)
``````

However, this is so computationally expensive way to match `SpatialLines` objects, especially if the number of objects is around 1000.

Is there any clever way to do this job?

P.S. Moreover, the above described approach does not guarantee all possible rotations and flips.

P.S2. If line #1 is zoomed out with respect to line #2, the the difference between line #1 and #2 must still be equal to 0.

UPDATE:

Any truly general-purpose effective method will standardize the representations of the shapes so that they will not change upon rotation, translation, reflection, or trivial changes in internal representation.

One way to do this is to list each connected shape as an alternating sequence of edge lengths and (signed) angles, starting from one end. (The shape should be "clean" in the sense of having no zero-length edges or straight angles.) To make this invariant under reflection, negate all the angles if the first nonzero one is negative.

(Because any connected polyline of n vertices will have n-1 edges separated by n-2 angles, I have found it convenient in the `R` code below to use a data structure consisting of two arrays, one for the edge lengths `\$lengths` and the other for the angles, `\$angles`. A line segment will have no angles at all, so it's important to handle zero-length arrays in such a data structure.)

Such representations can be ordered lexicographically. Some allowance should be made for floating-point errors accumulated during the standardization process. An elegant procedure would estimate those errors as a function of the original coordinates. In the solution below, a simpler method is used in which two lengths are considered equal when they differ by a very small amount on a relative basis. Angles may differ only by a very small amount on an absolute basis.

To make them invariant under reversal of the underlying orientation, choose the lexicographically earliest representation between that of the polyline and its reversal.

To handle multi-part polylines, arrange their components in lexicographic order.

To find the equivalence classes under Euclidean transformations, then,

• Create standardized representations of the shapes.

• Perform a lexicographic sort of the standardized representations.

• Make a pass through the sorted order to identify sequences of equal representations.

The computation time is proportional to O(n*log(n)*N) where n is the number of features and N is the largest number of vertices in any feature. This is efficient.

It's probably worth mentioning in passing that a preliminary grouping based on easily-calculated invariant geometric properties, such as polyline length, center, and moments about that center, often can be applied to streamline the entire process. One only needs to find subgroups of congruent features within each such preliminary group. The full method given here would be needed for shapes that otherwise would be so remarkably similar that such simple invariants still would not distinguish them. Simple features constructed from raster data might have such characteristics, for instance. However, since the solution given here is so efficient anyway, that if one is going to go to the effort of implementing it, it might work just fine all by itself.

Example

The left hand figure shows five polylines plus 15 more that were obtained from those via random translation, rotation, reflection, and reversal of internal orientation (which is not visible). The right hand figure colors them according to their Euclidean equivalence class: all figures of a common color are congruent; different colors are not congruent.

`R` code follows. When the inputs were updated to 500 shapes, 500 extra (congruent) shapes, with a mean of 100 vertices per shape, the execution time on this machine was 3 seconds.

This code is incomplete: because `R` does not have a native lexicographic sort, and I did not feel like coding one from scratch, I simply perform the sorting on the first coordinate of each standardized shape. That will be fine for the random shapes created here, but for production work a full lexicographic sort should be implemented. The function `order.shape` would be the only one affected by this change. Its input is a list of standardized shape `s` and its output is the sequence of indexes into `s` that would sort it.

``````#
# Create random shapes.
#
n.shapes <- 5      # Unique shapes, up to congruence
n.shapes.new <- 15 # Additional congruent shapes to generate
p.mean <- 5        # Expected number of vertices per shape
set.seed(17)       # Create a reproducible starting point
shape.random <- function(n) matrix(rnorm(2*n), nrow=2, ncol=n)
shapes <- lapply(2+rpois(n.shapes, p.mean-2), shape.random)
#
# Randomly move them around.
#
move.random <- function(xy) {
a <- runif(1, 0, 2*pi)
reflection <- sign(runif(1, -1, 1))
translation <- runif(2, -8, 8)
m <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2) %*%
matrix(c(reflection, 0, 0, 1), 2, 2)
m <- m %*% xy + translation
if (runif(1, -1, 0) < 0) m <- m[ ,dim(m)[2]:1]
return (m)
}
i <- sample(length(shapes), n.shapes.new, replace=TRUE)
shapes <- c(shapes, lapply(i, function(j) move.random(shapes[[j]])))
#
# Plot the shapes.
#
range.shapes <- c(min(sapply(shapes, min)), max(sapply(shapes, max)))
palette(gray.colors(length(shapes)))
par(mfrow=c(1,2))
plot(range.shapes, range.shapes, type="n",asp=1, bty="n", xlab="", ylab="")
invisible(lapply(1:length(shapes), function(i) lines(t(shapes[[i]]), col=i, lwd=2)))
#
# Standardize the shape description.
#
standardize <- function(xy) {
n <- dim(xy)[2]
vectors <- xy[ ,-1, drop=FALSE] - xy[ ,-n, drop=FALSE]
lengths <- sqrt(colSums(vectors^2))
if (which.min(lengths - rev(lengths))*2 < n) {
lengths <- rev(lengths)
vectors <- vectors[, (n-1):1]
}
if (n > 2) {
vectors <- vectors / rbind(lengths, lengths)
perps <- rbind(-vectors[2, ], vectors[1, ])
angles <- sapply(1:(n-2), function(i) {
cosine <- sum(vectors[, i+1] * vectors[, i])
sine <- sum(perps[, i+1] * vectors[, i])
atan2(sine, cosine)
})
i <- min(which(angles != 0))
angles <- sign(angles[i]) * angles
} else angles <- numeric(0)
list(lengths=lengths, angles=angles)
}
shapes.std <- lapply(shapes, standardize)
#
# Sort lexicographically.  (Not implemented: see the text.)
#
order.shape <- function(s) {
order(sapply(s, function(s) s\$lengths[1]))
}
i <- order.shape(shapes.std)
#
# Group.
#
equal.shape <- function(s.0, s.1) {
same.length <- function(a,b) abs(a-b) <= (a+b) * 1e-8
same.angle <- function(a,b) min(abs(a-b), abs(a-b)-2*pi) < 1e-11
r <- function(u) {
a <- u\$angles
if (length(a) > 0) {
a <- rev(u\$angles)
i <- min(which(a != 0))
a <- sign(a[i]) * a
}
list(lengths=rev(u\$lengths), angles=a)
}
e <- function(u, v) {
if (length(u\$lengths) != length(v\$lengths)) return (FALSE)
all(mapply(same.length, u\$lengths, v\$lengths)) &&
all(mapply(same.angle, u\$angles, v\$angles))
}
e(s.0, s.1) || e(r(s.0), s.1)
}
g <- rep(1, length(shapes.std))
for (j in 2:length(i)) {
i.0 <- i[j-1]
i.1 <- i[j]
if (equal.shape(shapes.std[[i.0]], shapes.std[[i.1]]))
g[j] <- g[j-1] else g[j] <- g[j-1]+1
}
palette(rainbow(max(g)))
plot(range.shapes, range.shapes, type="n",asp=1, bty="n", xlab="", ylab="")
invisible(lapply(1:length(i), function(j) lines(t(shapes[[i[j]]]), col=g[j], lwd=2)))
``````
• When one includes arbitrary dilations (or "isotheties") in the group of transformations, the equivalence classes are the congruence classes of affine geometry. This complication is easily handled: standardize all polylines to have total unit length, for instance. Jan 19 '15 at 22:28
• Thanks very much. Just one question: Should shapes be represented as SpatialLines or SpatialPolygons? Jan 19 '15 at 22:38
• Polygons create another complication: their boundaries have no definite endpoints. There are many ways to handle that, such as standardizing the representation to begin at (say) the vertex that sorts first in x-y lexicographic order and proceeding in a counterclockwise direction around the polygon. (A topologically "clean" connected polygon will have only one such vertex.) Whether a shape is considered a polygon or polyline depends on what kind of feature it represents: there is no intrinsic way to say of any closed list of points whether it is intended to be a polyline or polygon. Jan 19 '15 at 22:40
• Sorry for a simple question, but I should ask it to understand your example. Your object shapes.std has both \$lengths and \$angles. If, however, I run this code on my xy data (e.g. [1,] 3093.5 -2987.8 [2,] 3072.7 -2991.0 etc), it does not estimate angles, neither draws shapes. If I run plot(shapes[[1]]), then I can clearly see my polyline. So, how should I save polylines in R to be able to test your code on my data? Jan 19 '15 at 22:59
• I started with the same data structure you did: an array of (x,y) coordinates. My arrays put those coordinates into columns (as if you had used `rbind(x,y)` instead of `cbind(x,y)`). That's all you need: the `sp` library is not used. If you would like to follow what is done in detail, I suggest you start out with, say, `n.shapes <- 2`, `n.shapes.new <- 3`, and `p.mean <- 1`. Then `shapes`, `shapes.std`, etc. are all small enough to be easily inspected. The elegant--and "right" way--to deal with all this would be to create a class of standardized feature representations. Jan 19 '15 at 23:03

You're asking for a lot with arbitrary rotation and dilation! Not sure how useful Hausdorff distance would be there, but check it out. My approach would be reducing the number of cases to check via cheap data. For example, you could skip expensive comparisons if the length of the two linestrings is not an integer ratio (assuming integer/graduated scaling). You could similarly check if the bounding box area or their convex hull areas are in a nice ratio. I'm sure there are plenty of cheap checks you could do against the centroid, like distances or angles from start/end.

Only then, if you detect scaling, undo it and do the really expensive checks.

Clarification: I don't know the packages you are using. By integer ratio I meant you should divide both distances, check if the result is an integer, if not, invert that value (could be you picked the wrong order) and recheck. If you get an integer or close enough, you can deduce that there was perhaps scaling going on. Or it could just be two different shapes.

As for the bounding box, you probably got opposite points of the rectangle that represents it, so getting the area out of them is simple arithmetic. The principle behind the ratio comparison is the same, just that the result would be squared. Don't bother with convex hulls if you can't get them out of that R package nicely, it was just an idea (likely not cheap enough anyway).

• Thanks a lot. Could you please explain how to detect if the length of the two linestrings is not an integer ratio? Also, I appreciate a lot if you can give an example of checking "if the bounding box area or convex hull areas are in a nice ratio" Jan 19 '15 at 20:35
• For instance, if I extract spatial bounding box from spatial data, then I just receive two points: spl <- sp::SpatialLines(list(Lines(Line(xy.sp), ID=i))) b <- bbox(spl) Jan 19 '15 at 20:50
• Extended the main post. Jan 19 '15 at 21:24
• "If you get an integer or close enough, you can deduce that there was perhaps scaling going on." Couldn't a user have applied a scale of 1.4 or so? Jan 19 '15 at 21:35
• Sure, but my assumption was made clear, especially with later edits. I was imagining webmap style zooming, where one is nicely limited. Jan 19 '15 at 21:54

A good method to compare these polylines would be to rely on a representation as a sequence of (distances, turn angles) at each vertice: For a line composed of points `P1, P2, ..., PN`, such sequence would be:

(distance(P1P2), angle(P1,P2,P3), distance(P2P3),... ,angle(P(N-2),P(N-1),PN), distance(P(N-1)PN)).

According to your requirements, two lines are equals if and only if their corresponding sequences are the same (modulo the order and angle direction). Comparing number sequences is trivial.

By computing each polyline sequence only once and, as suggested by lynxlynxlynx, testing the sequence similarity only for polyline having same trivial characteristics (length, vertices number...), the computation should be really fast!

• This is the right idea. For it actually to work, though, many details need to be addressed, such as coping with reflections, internal orientation, possibility of multiple connected components, and floating point roundoff error. They are discussed in the solution I provided. Jan 19 '15 at 22:29
• Yes, I only described the main idea. Your answer is remarkably more complete (like often :-) Jan 20 '15 at 8:41