Trouver un rectangle minimum pour des points donnés?

72

Comme vous le voyez sur la figure, la question est la suivante:

Comment trouver le rectangle minimum de surface (MAR) ajusté sur les points donnés?

et une question complémentaire est:

Existe-t-il une solution analytique au problème?

(Un développement de la question consistera à adapter une boîte (3D) à un groupe de points dans un nuage de points 3D.)

Dans un premier temps, je propose de trouver la coque convexe pour les points qui reforment le problème (en supprimant ces points ne sont pas impliqués dans la solution): lier un MAR à un polygone. La méthode requise fournira X ( centre du rectangle ), D ( deux dimensions ) et A ( angle ).


Ma proposition de solution:

  • Trouver le centroïde du polygone (voir Recherche du centre de la géométrie de l'objet? )
  • [S] Ajuster un rectangle ajusté simple, c'est-à-dire parallèle aux axes X et Y
    • vous pouvez utiliser la minmaxfonction pour X et Y des points donnés (par exemple, les sommets d'un polygone)
  • Stocker la surface du rectangle ajusté
  • Faire pivoter le polygone autour du centre de gravité, par exemple, d'un degré
  • Répéter de [S] jusqu'à une rotation complète
  • Signaler l'angle de l'aire minimale comme résultat

Cela me semble prometteur, mais les problèmes suivants existent:

  • choisir une bonne résolution pour le changement d’angle pourrait être difficile,
  • le coût de calcul est élevé,
  • la solution n'est pas analytique mais expérimentale.

entrez la description de l'image ici

Développeur
la source

Réponses:

45

Oui, il existe une solution analytique à ce problème. L'algorithme que vous recherchez est connu dans la généralisation des polygones sous le nom de "plus petit rectangle environnant".

L'algorithme que vous décrivez est correct, mais pour résoudre les problèmes énumérés, vous pouvez utiliser le fait que l'orientation du MAR est identique à celle de l'un des bords de la coque convexe du nuage de points . Il suffit donc de tester les orientations des bords de la coque convexe. Vous devriez:

  • Calcule la coque convexe du nuage.
  • Pour chaque bord de la coque convexe:
    • calculer l'orientation des arêtes (avec arctan),
    • faire pivoter la coque convexe en utilisant cette orientation afin de calculer facilement la surface du rectangle de délimitation avec les min / max de x / y de la coque convexe tournée,
    • Enregistrer l'orientation correspondant à la surface minimale trouvée,
  • Renvoie le rectangle correspondant à la surface minimale trouvée.

Un exemple d'implémentation en Java est disponible ici .

En 3D, il en va de même, sauf que:

  • La coque convexe sera un volume,
  • Les orientations testées seront les orientations (en 3D) des faces de la coque convexe.

Bonne chance!

Julien
la source
11
+1 Très belle réponse! Je tiens à souligner que la rotation du nuage n’est pas nécessaire. Premièrement - vous vouliez probablement dire cela - seuls les sommets de la coque doivent être pris en compte. Deuxièmement, au lieu de tourner, représentez le côté actuel sous forme d'une paire de vecteurs unitaires orthogonaux. Le fait de prendre leurs produits de points avec les coordonnées de sommet de coque (ce qui pourrait être fait en une seule opération de matrice) donne les coordonnées pivotées: aucune trigonométrie n'est nécessaire, rapide et parfaitement précise.
whuber
2
Merci pour les liens. En effet, la rotation uniquement pour le nombre d'arêtes rend la méthode proposée très efficace. Je pourrais trouver le papier prouve que. Bien que j’ai marqué cela comme la réponse à la loyauté à la première bonne réponse (je ne peux pas choisir deux / plus de bonnes réponses :(), je voudrais recommander fortement de prendre en compte la réponse complète de whuber ci-dessous. L’efficacité de la méthode donnée ici (en évitant les rotations!) Est incroyable, et la procédure n’est constituée que de quelques lignes de code, elle est facilement traduisible en Python :)
Développeur
Pouvez-vous s'il vous plaît mettre à jour le lien d'implémentation Java?
Myra
oui c'est fait!
Julien
1
Notez que l'extension en 3D est un peu plus compliquée que cela. Chaque face de la coque convexe 3D définit une orientation possible d' une face du cadre de sélection, mais pas l'orientation des faces qui lui sont perpendiculaires. Le problème de la rotation du rectangle dans ce plan devient le problème 2D rectangle de délimitation minimal dans le plan de cette face. Pour chaque bord de la coque convexe du nuage projeté sur un plan donné, vous pouvez dessiner un cadre de sélection qui vous donnera un volume différent en 3D.
Sera
40

Pour compléter l'excellente solution de @ julien, voici une implémentation fonctionnelle dans R, qui pourrait servir de pseudocode pour guider toute implémentation spécifique au SIG (ou être appliquée directement dans R, bien sûr). L'entrée est un tableau de coordonnées de points. La sortie (la valeur de mbr) est un tableau des sommets du rectangle de délimitation minimal (le premier étant répété pour le fermer). Notez l'absence complète de tout calcul trigonométrique.

MBR <- function(p) {
  # Analyze the convex hull edges     
  a <- chull(p)                                   # Indexes of extremal points
  a <- c(a, a[1])                                 # Close the loop
  e <- p[a[-1],] - p[a[-length(a)], ]             # Edge directions
  norms <- sqrt(rowSums(e^2))                     # Edge lengths
  v <- e / norms                                  # Unit edge directions
  w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

  # Find the MBR
  vertices <- p[a, ]                              # Convex hull vertices
  x <- apply(vertices %*% t(v), 2, range)         # Extremes along edges
  y <- apply(vertices %*% t(w), 2, range)         # Extremes normal to edges
  areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
  k <- which.min(areas)                           # Index of the best edge (smallest area)

  # Form a rectangle from the extremes of the best edge
  cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}

Voici un exemple d'utilisation:

# Create sample data
set.seed(23)
p <- matrix(rnorm(20*2), ncol=2)                 # Random (normally distributed) points
mbr <- MBR(points)

# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, range) # Plotting limits
plot(p[(function(x) c(x, x[1]))(chull(p)), ], 
     type="l", asp=1, bty="n", xaxt="n", yaxt="n",
     col="Gray", pch=20, 
     xlab="", ylab="",
     xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=3)                         # The MBR
points(points, pch=19)                                # The points

MBR

Le timing est limité par la vitesse de l'algorithme de coque convexe, car le nombre de sommets dans la coque est presque toujours beaucoup plus petit que le total. La plupart des algorithmes de coque convexe sont asymptotiquement O (n * log (n)) pour n points: vous pouvez calculer presque aussi vite que vous pouvez lire les coordonnées.

whuber
la source
+1 Quelle solution incroyable! Une telle idée ne vient qu'après de longues expériences. À partir de maintenant, je serai curieux d’optimiser mes codes existants en s’inspirant de cette excellente réponse.
Développeur
Je souhaite que je pourrais upvote cela deux fois. J'apprends R et vos réponses sont une source d'inspiration continue.
John Powell
1
@retrovius Le rectangle de délimitation d'un ensemble de points (pivotés) est déterminé par quatre nombres: la plus petite coordonnée x, la plus grande coordonnée x, la plus petite coordonnée y et la plus grande coordonnée y. C'est ce que les "extrêmes le long des bords" se réfère.
whuber
1
@retrovius L'origine ne joue aucun rôle dans ces calculs, car tout est basé sur des différences de coordonnées, sauf à la fin, où le meilleur rectangle calculé en coordonnées pivotées est simplement retourné. Bien qu'il soit judicieux d'utiliser un système de coordonnées dans lequel l'origine est proche des points (pour minimiser la perte de précision en virgule flottante), l'origine n'est pas pertinente autrement.
whuber
1
@Retrovius Vous pouvez interpréter cela en termes d'une propriété de rotations: la matrice d'une rotation est orthogonale. Ainsi, un type de ressource serait une étude de l’algèbre linéaire (en général) ou de la géométrie euclidienne analytique (en particulier). Cependant, j’ai trouvé que le moyen le plus simple de gérer les rotations (et les translations et les redimensionnements) dans le plan est de considérer les points comme des nombres complexes: les rotations sont simplement effectuées en multipliant les valeurs par des nombres de longueurs unitaires.
whuber
8

Je viens juste de le mettre en œuvre moi-même et de poster ma réponse sur StackOverflow , mais je me suis dit que je laisserais tomber ma version ici pour que les autres le voient:

import numpy as np
from scipy.spatial import ConvexHull

def minimum_bounding_rectangle(points):
    """
    Find the smallest bounding rectangle for a set of points.
    Returns a set of points representing the corners of the bounding box.

    :param points: an nx2 matrix of coordinates
    :rval: an nx2 matrix of coordinates
    """
    from scipy.ndimage.interpolation import rotate
    pi2 = np.pi/2.

    # get the convex hull for the points
    hull_points = points[ConvexHull(points).vertices]

    # calculate edge angles
    edges = np.zeros((len(hull_points)-1, 2))
    edges = hull_points[1:] - hull_points[:-1]

    angles = np.zeros((len(edges)))
    angles = np.arctan2(edges[:, 1], edges[:, 0])

    angles = np.abs(np.mod(angles, pi2))
    angles = np.unique(angles)

    # find rotation matrices
    # XXX both work
    rotations = np.vstack([
        np.cos(angles),
        np.cos(angles-pi2),
        np.cos(angles+pi2),
        np.cos(angles)]).T
#     rotations = np.vstack([
#         np.cos(angles),
#         -np.sin(angles),
#         np.sin(angles),
#         np.cos(angles)]).T
    rotations = rotations.reshape((-1, 2, 2))

    # apply rotations to the hull
    rot_points = np.dot(rotations, hull_points.T)

    # find the bounding points
    min_x = np.nanmin(rot_points[:, 0], axis=1)
    max_x = np.nanmax(rot_points[:, 0], axis=1)
    min_y = np.nanmin(rot_points[:, 1], axis=1)
    max_y = np.nanmax(rot_points[:, 1], axis=1)

    # find the box with the best area
    areas = (max_x - min_x) * (max_y - min_y)
    best_idx = np.argmin(areas)

    # return the best box
    x1 = max_x[best_idx]
    x2 = min_x[best_idx]
    y1 = max_y[best_idx]
    y2 = min_y[best_idx]
    r = rotations[best_idx]

    rval = np.zeros((4, 2))
    rval[0] = np.dot([x1, y2], r)
    rval[1] = np.dot([x2, y2], r)
    rval[2] = np.dot([x2, y1], r)
    rval[3] = np.dot([x1, y1], r)

    return rval

Voici quatre exemples différents de celui-ci en action. Pour chaque exemple, j'ai généré 4 points aléatoires et trouvé la boîte englobante.

entrez la description de l'image ici

C'est relativement rapide aussi pour ces échantillons sur 4 points:

>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop
JesseBuesking
la source
Bonjour JesseBuesking, pouvez-vous générer des rectangles avec des angles de 90 degrés? Votre code fonctionne bien pour obtenir des parallélogrammes, mais des angles de 90 degrés sont nécessaires dans mon cas d'utilisation spécifique. Pourriez-vous recommander comment votre code peut être modifié pour atteindre cela? Merci!
Nader Alexan
@NaderAlexan Si vous vous demandez s'il peut gérer les carrés, alors oui, c'est le cas! Je viens de l'essayer sur une unité carrée points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]]), et le résultat array([[1.00000000e+00, 6.12323400e-17], [0.00000000e+00, 0.00000000e+00], [6.12323400e-17, 1.00000000e+00], [1.00000000e+00, 1.00000000e+00]])obtenu correspond à l'unité carrée (y compris certaines erreurs d'arrondi en virgule flottante). Note: un carré est juste un rectangle avec des côtés égaux, donc je suppose que s'il peut gérer un carré, il se généralise à tous les rectangles.
JesseBuesking
Merci pour votre réponse. Oui, cela fonctionne très bien, mais j'essaie de le forcer à toujours produire un rectangle (4 côtés avec des angles de 90 degrés pour chaque côté) par rapport à tout autre polygone à 4 côtés, bien que dans certains cas cela produise un rectangle, cela ne semble pas Pour être une contrainte constante, savez-vous comment modifier le code pour ajouter cette contrainte? Merci!
Nader Alexan
Peut-être que gis.stackexchange.com/a/22934/48041 pourrait vous guider vers une solution, étant donné que leur réponse semble avoir cette contrainte? Une fois que vous avez trouvé une solution, vous devriez y contribuer, car je suis sûr que d’autres le trouveront utile. Bonne chance!
JesseBuesking
7

Il existe un outil dans Whitebox GAT ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ) appelé Minimum Bounding Box pour résoudre ce problème précis. Il y a aussi un outil de coque convexe minimum là aussi. Plusieurs outils de la boîte à outils Patch Shape, par exemple l'orientation et l'allongement du patch, sont basés sur la recherche du cadre de sélection minimal.

entrez la description de l'image ici


la source
4

Je suis tombé sur ce fil en cherchant une solution Python pour un rectangle de délimitation à aire minimale.

Voici ma mise en œuvre , pour laquelle les résultats ont été vérifiés avec Matlab.

Le code de test est inclus pour les polygones simples, et je l’utilise pour trouver le cadre de délimitation minimal 2D et les directions des axes pour un PointCloud 3D.

David
la source
Votre réponse a-t-elle été supprimée?
Paul Richter
@ PaulRichter apparemment. La source était ici github.com/dbworth/minimum-area-bounding-rectangle bien
sehe
3

Merci @ Whuber réponse. C'est une excellente solution, mais lente pour les gros nuages ​​de points. J'ai trouvé que la convhullnfonction dans le package R geometryest beaucoup plus rapide (138 s contre 0,03 s pour 200000 points). J'ai collé mes codes ici pour quiconque est intéressant pour une solution plus rapide.

library(alphahull)                                  # Exposes ashape()
MBR <- function(points) {
    # Analyze the convex hull edges                       
    a <- ashape(points, alpha=1000)                 # One way to get a convex hull...
    e <- a$edges[, 5:6] - a$edges[, 3:4]            # Edge directions
    norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths
    v <- diag(1/norms) %*% e                        # Unit edge directions
    w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

    # Find the MBR
    vertices <- (points) [a$alpha.extremes, 1:2]    # Convex hull vertices
    minmax <- function(x) c(min(x), max(x))         # Computes min and max
    x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
    y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
    areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
    k <- which.min(areas)                           # Index of the best edge (smallest area)

    # Form a rectangle from the extremes of the best edge
    cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,])
}

MBR2 <- function(points) {
    tryCatch({
        a2 <- geometry::convhulln(points, options = 'FA')

        e <- points[a2$hull[,2],] - points[a2$hull[,1],]            # Edge directions
        norms <- apply(e, 1, function(x) sqrt(x %*% x)) # Edge lengths

        v <- diag(1/norms) %*% as.matrix(e)                        # Unit edge directions


        w <- cbind(-v[,2], v[,1])                       # Normal directions to the edges

        # Find the MBR
        vertices <- as.matrix((points) [a2$hull, 1:2])    # Convex hull vertices
        minmax <- function(x) c(min(x), max(x))         # Computes min and max
        x <- apply(vertices %*% t(v), 2, minmax)        # Extremes along edges
        y <- apply(vertices %*% t(w), 2, minmax)        # Extremes normal to edges
        areas <- (y[1,]-y[2,])*(x[1,]-x[2,])            # Areas
        k <- which.min(areas)                           # Index of the best edge (smallest area)

        # Form a rectangle from the extremes of the best edge
        as.data.frame(cbind(x[c(1,2,2,1,1),k], y[c(1,1,2,2,1),k]) %*% rbind(v[k,], w[k,]))
    }, error = function(e) {
        assign('points', points, .GlobalEnv)
        stop(e)  
    })
}


# Create sample data
#set.seed(23)
points <- matrix(rnorm(200000*2), ncol=2)                 # Random (normally distributed) points
system.time(mbr <- MBR(points))
system.time(mmbr2 <- MBR2(points))


# Plot the hull, the MBR, and the points
limits <- apply(mbr, 2, function(x) c(min(x),max(x))) # Plotting limits
plot(ashape(points, alpha=1000), col="Gray", pch=20, 
     xlim=limits[,1], ylim=limits[,2])                # The hull
lines(mbr, col="Blue", lwd=10)                         # The MBR
lines(mbr2, col="red", lwd=3)                         # The MBR2
points(points, pch=19)   

Deux méthodes obtiennent la même réponse (exemple pour 2000 points):

entrez la description de l'image ici

Bangyou
la source
Est-il possible d'étendre cette implémentation à un espace 3D (c'est-à-dire à trouver une boîte de volume minimum incluant tous les points donnés dans un espace 3D)?
Sasha
0

Je recommande simplement la fonction intégrée de OpenCV minAreaRect, qui trouve un rectangle pivoté de la surface minimale entourant le jeu de points 2D en entrée. Pour voir comment utiliser cette fonction, on peut se référer à ce tutoriel .

piège
la source