Znalezienie prostokąta o minimalnej powierzchni dla danych punktów?

Jak widać na rysunku, pytanie brzmi:

Jak znaleźć prostokąt o minimalnej powierzchni (MAR) dopasowany do danych punktów?

a pytanie pomocnicze brzmi:

Czy istnieje jakieś analityczne rozwiązanie problemu?

(Opracowanie pytania będzie polegało na dopasowaniu ramki (3D) do klastra punktów w chmurze punktów 3D).

Jako pierwszy etap proponuję znaleźć wypukły kadłub dla punktów, które reformują problem (poprzez usunięcie tych punktów nie są zaangażowane w rozwiązanie) w celu: dopasowania MAR do wielokąta. Wymagana metoda zapewni X ( środek prostokąta ), D ( dwa wymiary ) i A ( kąt ).

Moja propozycja rozwiązania:

• Znajdź środek ciężkości wielokąta (patrz Znajdowanie środka geometrii obiektu? )
• [S] Dopasuj prosty dopasowany prostokąt, tzn. Równolegle do osi X i Y
  • możesz użyć minmaxfunkcji dla X i Y podanych punktów (np. wierzchołków wielokąta)
• Przechowuj obszar dopasowanego prostokąta
• Obróć wielokąt wokół środka ciężkości o np. 1 stopień
• Powtarzaj od [S] aż do pełnego obrotu
• Podaj jako wynik kąt minimalnej powierzchni

Wydaje mi się to obiecujące, jednak istnieją następujące problemy:

• wybór dobrej rozdzielczości dla zmiany kąta może być trudny,
• koszt obliczeń jest wysoki,
• rozwiązanie nie jest analityczne, ale eksperymentalne.

Tak, istnieje analityczne rozwiązanie tego problemu. Algorytm, którego szukasz, jest znany w uogólnieniu wielokąta jako „najmniejszy otaczający prostokąt”.

Algorytm, który opisujesz, jest w porządku, ale aby rozwiązać wymienione problemy, możesz wykorzystać fakt, że orientacja MAR jest taka sama jak jedna z krawędzi wypukłego kadłuba chmury punktów . Wystarczy więc przetestować orientację wypukłych krawędzi kadłuba. Powinieneś:

• Oblicz wypukły kadłub chmury.
• Dla każdej krawędzi wypukłego kadłuba:
  • oblicz orientację krawędzi (z arctan),
  • obrócić wypukły kadłub za pomocą tej orientacji, aby łatwo obliczyć obwiednię prostokąta z min / max x / y obróconego wypukłego kadłuba,
  • Przechowuj orientację odpowiadającą minimalnemu znalezionemu obszarowi,
• Zwraca prostokąt odpowiadający znalezionemu minimalnemu obszarowi.

Przykład implementacji w Javie jest tam dostępny .

W 3D to samo dotyczy, z wyjątkiem:

• Wypukły kadłub będzie objętością,
• Testowane orientacje będą orientacjami (w 3D) wypukłych powierzchni kadłuba.

Powodzenia!

Aby uzupełnić świetne rozwiązanie @ julien, oto działająca implementacja w R, która może służyć jako pseudokod, który poprowadzi każdą implementację specyficzną dla GIS (lub Roczywiście można ją zastosować bezpośrednio ). Dane wejściowe to tablica współrzędnych punktowych. Dane wyjściowe (wartość mbr) to tablica wierzchołków minimalnego prostokąta ograniczającego (z pierwszym powtórzeniem w celu zamknięcia). Zwróć uwagę na całkowity brak jakichkolwiek obliczeń trygonometrycznych.

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,])
}

Oto przykład jego użycia:

# 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

Czas jest ograniczony przez szybkość algorytmu wypukłego kadłuba, ponieważ liczba wierzchołków kadłuba jest prawie zawsze znacznie mniejsza niż suma. Większość wypukłych algorytmów kadłuba ma asymptotycznie O (n * log (n)) dla n punktów: można obliczyć prawie tak szybko, jak można odczytać współrzędne.

Właśnie to zaimplementowałem i opublikowałem swoją odpowiedź na StackOverflow , ale pomyślałem, że upuszczę moją wersję tutaj, aby inni mogli ją zobaczyć:

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

Oto cztery różne przykłady tego działania. Dla każdego przykładu wygenerowałem 4 losowe punkty i znalazłem obwiednię.

Dla tych próbek jest też stosunkowo szybki w 4 punktach:

>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop

W Whitebox GAT ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ) znajduje się narzędzie o nazwie Minimalna ramka ograniczająca do rozwiązania tego konkretnego problemu. Jest tam również narzędzie z minimalnym wypukłym kadłubem. Kilka narzędzi w przyborniku kształtu łatki, np. Orientacja i wydłużenie łatki, opiera się na znalezieniu minimalnego obwiedni.

Natknąłem się na ten wątek, szukając rozwiązania w języku Python dla prostokąta ograniczającego obszar minimalny.

Oto moja implementacja , dla której wyniki zostały zweryfikowane za pomocą Matlaba.

Dołączono kod testowy dla prostych wielokątów i używam go do znalezienia minimalnej ramki granicznej 2D i kierunków osi dla 3D PointCloud.

Dzięki @ odpowiedź Whubera. To świetne rozwiązanie, ale wolne dla dużych chmur punktów. Odkryłem, że convhullnfunkcja w pakiecie R geometryjest znacznie szybsza (138 s vs 0,03 s za 200000 punktów). Wkleiłem tutaj moje kody, aby każdy był ciekawy dla szybszego rozwiązania.

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)   

Dwie metody uzyskują tę samą odpowiedź (przykład na 2000 punktów):

Po prostu polecam wbudowaną funkcję OpenCV minAreaRect, która znajduje obrócony prostokąt o minimalnym obszarze otaczającym wejściowy zestaw punktów 2D. Aby zobaczyć, jak korzystać z tej funkcji, zapoznaj się z tym samouczkiem .