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

71

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.

wprowadź opis zdjęcia tutaj

Deweloper
źródło

Odpowiedzi:

45

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!

Julien
źródło
11
+1 Bardzo ładna odpowiedź! Chciałbym zauważyć, że faktyczny obrót chmury jest niepotrzebny. Po pierwsze - prawdopodobnie miałeś na myśli to - należy brać pod uwagę tylko wierzchołki kadłuba. Po drugie, zamiast obracać, reprezentują stronę bieżącą jako parę wektorów jednostek ortogonalnych. Biorąc ich produkty kropkowe ze współrzędnymi wierzchołka kadłuba (co można zrobić jako operację pojedynczej matrycy), daje obrócone współrzędne: nie wymaga trygonometrii, jest szybka i idealnie dokładna.
whuber
2
Dzięki za linki. Rzeczywiście, obracanie tylko dla # krawędzi sprawia, że ​​proponowana metoda jest bardzo wydajna. Mogłem znaleźć artykuł, który to potwierdza. Chociaż zaznaczyłem to jako odpowiedź na lojalność w stosunku do pierwszej dobrej odpowiedzi (nie mogę wybrać dwóch / więcej świetnych odpowiedzi :() Chciałbym gorąco polecić, biorąc pod uwagę pełną odpowiedź Whubera poniżej. Wydajność tamtej metody (unikanie rotacji!) Wynosi niesamowite, a cała procedura to tylko kilka wierszy kodu. Dla mnie jest to łatwe do przetłumaczenia na Python :)
programista
Czy możesz zaktualizować link do implementacji Java?
Myra
tak, gotowe!
Julien
1
Pamiętaj, że rozszerzenie do 3D jest nieco bardziej skomplikowane. Każda powierzchnia wypukłego kadłuba 3D określa możliwą orientację jednej powierzchni obwiedni, ale nie orientację ścian prostopadłych do niej. Problem z obracaniem pudełka w tej płaszczyźnie staje się problemem prostokąta ograniczającego minimum 2D w płaszczyźnie tej powierzchni. Dla każdej krawędzi wypukłego kadłuba chmury rzutowanej na daną płaszczyznę możesz narysować obwiednię, która da ci inną objętość w 3D.
Czy
40

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

MBR

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.

Whuber
źródło
+1 Co za niesamowite rozwiązanie! Taki pomysł pojawia się dopiero po długich doświadczeniach. Odtąd będę ciekawy optymalizacji moich istniejących kodów, inspirowanych tą świetną odpowiedzią.
Deweloper
Chciałbym móc to zagłosować dwa razy. Uczę się R, a twoje odpowiedzi są ciągłym źródłem inspiracji.
John Powell,
1
@retrovius Granica prostokąta zbioru (obróconych) punktów jest wyznaczona przez cztery liczby: najmniejszą współrzędną x, największą współrzędną x, najmniejszą współrzędną y i największą współrzędną y. Właśnie do tego odnoszą się „skrajności wzdłuż krawędzi”.
whuber
1
@retrovius Początek nie odgrywa żadnej roli w tych obliczeniach, ponieważ wszystko opiera się na różnicach współrzędnych, z wyjątkiem końca, gdzie najlepszy prostokąt obliczony w obróconych współrzędnych jest po prostu obracany do tyłu. Chociaż rozsądnym pomysłem jest użycie układu współrzędnych, w którym punkt początkowy znajduje się blisko punktów (aby zminimalizować utratę precyzji zmiennoprzecinkowej), w przeciwnym razie początek nie jest istotny.
whuber
1
@Retrovius Można to interpretować w kategoriach właściwości obrotu: mianowicie macierz obrotu jest ortogonalna. Tak więc jednym rodzajem zasobów byłoby badanie algebry liniowej (ogólnie) lub analitycznej geometrii euklidesowej (konkretnie). Odkryłem jednak, że najłatwiejszym sposobem radzenia sobie z rotacjami (oraz translacjami i skalowaniem) w płaszczyźnie jest postrzeganie punktów jako liczb zespolonych: rotacje są po prostu wykonywane przez pomnożenie wartości przez liczby o długości jednostkowej.
whuber
8

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ę.

wprowadź opis zdjęcia tutaj

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
JesseBuesking
źródło
Cześć JesseBuesking, czy potrafisz generować prostokąty o narożnikach 90 stopni? Twój kod świetnie się sprawdza w uzyskiwaniu równoległoboków, ale w moim konkretnym przypadku użycia wymagane są narożniki 90 stopni. Czy możesz polecić, w jaki sposób można zmodyfikować kod, aby to osiągnąć? Dzięki!
Nader Alexan
@NaderAlexan Jeśli pytasz, czy poradzi sobie z kwadratami, to tak, na pewno tak! Właśnie wypróbowałem go na kwadracie jednostkowym points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]]), a wynik jest tym, 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]])który jest kwadratem jednostkowym (włączając pewne błędy zaokrąglania zmiennoprzecinkowego). Uwaga: kwadrat jest tylko prostokątem o równych bokach, więc zakładam, że może obsłużyć kwadrat, który uogólnia na wszystkie prostokąty.
JesseBuesking
Dziękuję za Twoją odpowiedź. Tak, działa świetnie, ale próbuję zmusić go, aby zawsze tworzył prostokąt (4 strony z kątami 90 stopni dla każdej strony) nad dowolnym innym 4-stronnym wielokątem, chociaż w niektórych przypadkach tworzy prostokąt, który nie wydaje się aby być stałym ograniczeniem, czy wiesz, jak zmodyfikować kod, aby dodać to ograniczenie? Dzięki!
Nader Alexan
Może gis.stackexchange.com/a/22934/48041 może poprowadzić cię w kierunku rozwiązania, biorąc pod uwagę, że ich odpowiedź wydaje się mieć takie ograniczenie? Gdy znajdziesz rozwiązanie, powinieneś je wnieść, ponieważ jestem pewien, że inni uznają je za przydatne. Powodzenia!
JesseBuesking
7

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.

wprowadź opis zdjęcia tutaj


źródło
4

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.

David
źródło
Czy twoja odpowiedź została usunięta?
Paul Richter,
@PaulRichter najwyraźniej. Źródłem tu github.com/dbworth/minimum-area-bounding-rectangle choć
sehe
3

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

wprowadź opis zdjęcia tutaj

Bangyou
źródło
Czy można rozszerzyć tę implementację na przestrzeń 3D (tj. Znaleźć pole minimalnej objętości, które obejmuje wszystkie podane punkty w przestrzeni 3D)?
Sasha,
0

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 .

pułapka
źródło