Menemukan minimum-area-persegi panjang untuk poin yang diberikan?

Seperti yang Anda lihat pada gambar, pertanyaannya adalah:

Bagaimana menemukan minimum-area-rectangle (MAR) terpasang pada poin yang diberikan?

dan pertanyaan pendukung adalah:

Apakah ada solusi analitis untuk masalah ini?

(Pengembangan pertanyaan adalah mencocokkan kotak (3D) ke sekelompok titik di awan titik 3D.)

Sebagai tahap pertama saya mengusulkan untuk menemukan cembung-hull untuk poin-poin yang mereformasi masalah (dengan menghilangkan titik-titik tersebut tidak terlibat dalam solusi) untuk: memasang MAR ke poligon. Metode yang diperlukan akan memberikan X ( pusat persegi panjang ), D ( dua dimensi ) dan A ( sudut ).

Proposal saya untuk solusi:

• Temukan pusat massa poligon (lihat Pusat penemuan geometri objek? )
• [S] Pas dengan persegi panjang sederhana, sejajar dengan sumbu X dan Y
  • Anda dapat menggunakan minmaxfungsi untuk X dan Y dari poin yang diberikan (mis., simpul poligon)
• Simpan area persegi panjang yang dipasang
• Putar poligon tentang centroid dengan misalnya 1 derajat
• Ulangi dari [S] hingga rotasi penuh selesai
• Laporkan sudut area minimum sebagai hasilnya

Bagi saya itu tampak menjanjikan, namun masalah berikut ada:

• memilih resolusi yang baik untuk perubahan sudut bisa jadi menantang,
• biaya perhitungannya tinggi,
• solusinya bukan analitis tetapi eksperimental.

Ya, ada solusi analitis untuk masalah ini. Algoritma yang Anda cari dikenal dalam generalisasi poligon sebagai "persegi panjang terkecil di sekitar".

Algoritma yang Anda gambarkan baik-baik saja tetapi untuk menyelesaikan masalah yang telah Anda daftarkan, Anda dapat menggunakan fakta bahwa orientasi MAR sama dengan salah satu dari salah satu ujung lambung cembung awan titik . Jadi, Anda hanya perlu menguji orientasi tepi lambung cembung. Anda harus:

• Hitung lambung cembung awan.
• Untuk setiap tepi lambung cembung:
  • menghitung orientasi tepi (dengan arctan),
  • putar cembung lambung menggunakan orientasi ini untuk menghitung dengan mudah daerah persegi panjang yang terikat dengan min / max x / y dari lambung cembung diputar,
  • Simpan orientasi sesuai dengan area minimum yang ditemukan,
• Kembalikan kotak yang sesuai dengan area minimum yang ditemukan.

Contoh implementasi di java tersedia di sana .

Dalam 3D, hal yang sama berlaku, kecuali:

• Lambung cembung akan menjadi volume,
• Orientasi yang diuji akan menjadi orientasi (dalam 3D) dari permukaan cembung cembung.

Semoga berhasil!

Untuk melengkapi solusi hebat @ julien, berikut ini adalah implementasi Ryang berfungsi, yang dapat berfungsi sebagai pseudocode untuk memandu implementasi spesifik GIS (atau diterapkan langsung R, tentu saja). Input adalah array koordinat titik. Output (nilai dari mbr) adalah array dari simpul persegi panjang batas minimum (dengan yang pertama diulang untuk menutupnya). Perhatikan tidak adanya perhitungan trigonometri.

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

Berikut ini adalah contoh penggunaannya:

# 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

Pengaturan waktu dibatasi oleh kecepatan algoritma cembung cembung, karena jumlah simpul dalam lambung hampir selalu jauh lebih sedikit daripada total. Sebagian besar algoritma cembung cangkang asimtotik O (n * log (n)) untuk n poin: Anda dapat menghitung hampir secepat Anda dapat membaca koordinat.

Saya baru saja mengimplementasikan ini sendiri dan memposting jawaban saya di StackOverflow , tetapi saya pikir saya akan menjatuhkan versi saya di sini agar orang lain dapat melihatnya:

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

Berikut ini empat contoh berbeda dalam aksi. Untuk setiap contoh, saya menghasilkan 4 poin acak dan menemukan kotak pembatas.

Ini juga relatif cepat untuk sampel-sampel ini pada 4 poin:

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

Ada alat di Whitebox GAT ( http://www.uoguelph.ca/~hydrogeo/Whitebox/ ) yang disebut Kotak Batas Minimum untuk menyelesaikan masalah yang pasti ini. Ada juga alat lambung cembung minimum di sana juga. Beberapa alat dalam kotak Bentuk Patch, misalnya orientasi dan perpanjangan patch, didasarkan pada menemukan kotak batas minimum.

Saya menemukan utas ini sambil mencari solusi Python untuk persegi panjang batas minimum.

Inilah implementasi saya , yang hasilnya diverifikasi dengan Matlab.

Kode uji disertakan untuk poligon sederhana, dan saya menggunakannya untuk menemukan kotak batas minimum 2D dan arah sumbu untuk 3D PointCloud.

Terima kasih @ whuber. Ini adalah solusi yang bagus, tetapi lambat untuk cloud big point. Saya menemukan convhullnfungsi dalam paket R geometryjauh lebih cepat (138 detik vs 0,03 detik untuk 200000 poin). Saya menempelkan kode saya di sini untuk siapa pun yang menarik untuk solusi yang lebih cepat.

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)   

Dua metode mendapatkan jawaban yang sama (contoh untuk 2000 poin):

Saya hanya merekomendasikan fungsi built-in OpenCV minAreaRect, yang menemukan persegi panjang yang diputar dari area minimum yang melampirkan input set poin 2D. Untuk melihat bagaimana menggunakan fungsi ini, orang dapat merujuk ke tutorial ini .