Comment mesurer la similitude des objets SpatialLines

J'ai créé deux SpatialLinesobjets R: .

Ces objets ont été créés de cette façon:

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

Maintenant, je veux en quelque sorte conclure qu'il s'agit de la même ligne tournée et inversée, et que la différence entre eux est égale à 0 (c'est-à-dire que la forme est égale).

Pour ce faire, on peut utiliser le maptoolspackage et faire pivoter la ligne # 1, par exemple:

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

Chaque ligne pivotée doit ensuite être vérifiée par rapport à la ligne n ° 2 à l'aide du rgeospackage, par exemple:

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

Cependant, c'est un moyen si coûteux en termes de calcul de faire correspondre des SpatialLinesobjets, surtout si le nombre d'objets est d'environ 1000.

Existe-t-il un moyen intelligent de faire ce travail?

PS De plus, l'approche décrite ci-dessus ne garantit pas toutes les rotations et retournements possibles.

P.S2. Si la ligne # 1 est dézoomée par rapport à la ligne # 2, la différence entre la ligne # 1 et # 2 doit toujours être égale à 0.

MISE À JOUR:

Toute méthode efficace vraiment polyvalente normalisera les représentations des formes afin qu'elles ne changent pas lors de la rotation, de la translation, de la réflexion ou de changements triviaux dans la représentation interne.

Une façon de procéder consiste à répertorier chaque forme connectée sous la forme d'une séquence alternée de longueurs de bord et d'angles (signés), en commençant par une extrémité. (La forme doit être "propre" dans le sens de ne pas avoir d'arêtes de longueur nulle ou d'angles droits.) Pour rendre cet invariant en réflexion, annulez tous les angles si le premier non nul est négatif.

(Parce que toute polyligne connectée de n sommets aura n -1 arêtes séparées par n -2 angles, j'ai trouvé pratique dans le Rcode ci-dessous d'utiliser une structure de données composée de deux tableaux, l'un pour les longueurs d'arête $lengthset l'autre pour le angles $angles. Un segment de ligne n'aura aucun angle, il est donc important de gérer des tableaux de longueur nulle dans une telle structure de données.)

Ces représentations peuvent être ordonnées lexicographiquement. Il convient de tenir compte des erreurs à virgule flottante accumulées au cours du processus de normalisation. Une procédure élégante estimerait ces erreurs en fonction des coordonnées d'origine. Dans la solution ci-dessous, une méthode plus simple est utilisée dans laquelle deux longueurs sont considérées comme égales lorsqu'elles diffèrent d'une très petite quantité sur une base relative. Les angles ne peuvent différer que d'une très petite quantité sur une base absolue.

Pour les rendre invariants sous inversion de l'orientation sous-jacente, choisissez la représentation lexicographiquement la plus ancienne entre celle de la polyligne et son inversion.

Pour gérer les polylignes en plusieurs parties, organisez leurs composants dans l'ordre lexicographique.

Pour trouver les classes d'équivalence sous les transformations euclidiennes, alors,

• Créez des représentations standardisées des formes.

• Effectuer une sorte lexicographique des représentations standardisées.

• Effectuez un passage dans l'ordre trié pour identifier des séquences de représentations égales.

Le temps de calcul est proportionnel à O (n * log (n) * N) où n est le nombre d'entités et N est le plus grand nombre de sommets d'une entité. C'est efficace.

Il vaut probablement la peine de mentionner en passant qu'un regroupement préliminaire basé sur des propriétés géométriques invariantes faciles à calculer, telles que la longueur de la polyligne, le centre et les moments autour de ce centre, peut souvent être appliqué pour rationaliser l'ensemble du processus. Il suffit de trouver des sous-groupes de caractéristiques congruentes dans chacun de ces groupes préliminaires. La méthode complète donnée ici serait nécessaire pour des formes qui autrement seraient si remarquablement similaires que de tels invariants simples ne les distingueraient toujours pas. Les entités simples construites à partir de données raster peuvent avoir de telles caractéristiques, par exemple. Cependant, étant donné que la solution donnée ici est si efficace de toute façon, que si l'on va faire l'effort de la mettre en œuvre, cela pourrait très bien fonctionner tout seul.

Exemple

La figure de gauche montre cinq polylignes plus 15 autres qui ont été obtenues à partir de celles-ci par translation aléatoire, rotation, réflexion et inversion de l'orientation interne (qui n'est pas visible). La figure de droite les colore selon leur classe d'équivalence euclidienne: toutes les figures d'une couleur commune sont congruentes; différentes couleurs ne sont pas congruentes.

Rle code suit. Lorsque les entrées ont été mises à jour en 500 formes, 500 formes supplémentaires (congruentes), avec une moyenne de 100 sommets par forme, le temps d'exécution sur cette machine était de 3 secondes.

Ce code est incomplet: parce qu'il Rn'a pas de tri lexicographique natif, et que je n'ai pas envie d'en coder un à partir de zéro, j'effectue simplement le tri sur la première coordonnée de chaque forme standardisée. Ce sera bien pour les formes aléatoires créées ici, mais pour le travail de production, un tri lexicographique complet devrait être implémenté. La fonction order.shapeserait la seule affectée par ce changement. Son entrée est une liste de formes normalisées set sa sortie est la séquence d'index dans slaquelle elle serait triée.

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

Vous en demandez beaucoup avec une rotation et une dilatation arbitraires! Je ne sais pas à quel point la distance Hausdorff serait utile, mais vérifiez-la. Mon approche consisterait à réduire le nombre de cas à vérifier via des données bon marché. Par exemple, vous pouvez ignorer les comparaisons coûteuses si la longueur des deux chaînes de lignes n'est pas un rapport entier (en supposant une mise à l'échelle entière / graduée ). Vous pouvez également vérifier si la zone du cadre de délimitation ou leurs zones de coque convexes sont dans un bon rapport. Je suis sûr qu'il y a beaucoup de contrôles bon marché que vous pouvez faire contre le centroïde, comme les distances ou les angles de début / fin.

Alors seulement, si vous détectez la mise à l'échelle, annulez-la et effectuez les vérifications vraiment coûteuses.

Clarification: je ne connais pas les packages que vous utilisez. Par rapport entier, je voulais dire que vous devriez diviser les deux distances, vérifier si le résultat est un entier, sinon, inverser cette valeur (vous pourriez avoir choisi le mauvais ordre) et revérifier. Si vous obtenez un entier ou suffisamment proche, vous pouvez en déduire qu'il y avait peut-être une mise à l'échelle en cours. Ou il pourrait simplement s'agir de deux formes différentes.

En ce qui concerne le cadre de sélection, vous avez probablement obtenu des points opposés du rectangle qui le représente, donc en extraire la zone est une simple arithmétique. Le principe derrière la comparaison des ratios est le même, juste que le résultat serait au carré. Ne vous embêtez pas avec des coques convexes si vous ne pouvez pas les retirer de ce package R, c'était juste une idée (probablement pas assez bon marché de toute façon).

Une bonne méthode pour comparer ces polylignes serait de s'appuyer sur une représentation comme une séquence de (distances, angles de virage) à chaque sommet: pour une ligne composée de points P1, P2, ..., PN, une telle séquence serait:

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

Selon vos besoins, deux lignes sont égales si et seulement si leurs séquences correspondantes sont identiques (modulez l'ordre et la direction de l'angle). La comparaison des séquences de nombres est triviale.

En calculant chaque séquence de polyligne une seule fois et, comme suggéré par lynxlynxlynx, en testant la similitude de séquence uniquement pour une polyligne ayant les mêmes caractéristiques triviales (longueur, nombre de sommets ...), le calcul devrait être vraiment rapide!