4

I'd like to calculate the "Oriented Bounding Box" (OBB) which is the enclosing rectangle aligned to the longest extent of a point set (which is not necessarily aligned with the axes of the co-ordinate system).

The solution could be similar to Finding minimum-area-rectangle for given points? but I'm looking for a robust solution in PostGIS SQL (and PL/Python if necessary).

image

Taras
  • 32,823
  • 4
  • 66
  • 137
Stefan
  • 811
  • 7
  • 12
  • 1
    That's a hard one.. I've done it for polygons in C# (or was the VB.net?) but not for points. You could start with a convex hull and find the longest axis I suppose. Do you have any other geometry tools available or just python/PostGIS? What part of the linked question isn't suitable? – Michael Stimson Feb 06 '15 at 00:36
  • I'm looking for a fast implementation in PL/pgSQL and tought to reuse ST_ConvexHull(). – Stefan Feb 06 '15 at 00:52
  • From the Convex Hull do some rotations and pick the smallest box to cover all the points. The hull (geometry) might have an extent property which is already the bounding box then it's a matter of recording the dimensions and angle, picking the smallest box and rotating the central line using the reverse matrix. – Michael Stimson Feb 06 '15 at 00:58
  • Centre of the envelope is ( XMax + XMin ) / 2, ( YMax + YMin ) / 2. Rotation with a matrix will allow rotation more than 90 deg - I wouldn't assume symmetry as your points may not be symmetric around the centre point. – Michael Stimson Feb 06 '15 at 01:18
  • Ok. But I still think there's potential for optimization, especially for avoiding to rotate 360 degrees. What about starting with ST_LongestLine (using the same geometry twice as input, as mentioned in linked question)? Then, only +/- 90 degrees need to be rotated? – Stefan Feb 07 '15 at 14:39

3 Answers3

2

This seems to be peanuts for JTS. Let's hope that the method is ported into GEOS so that you can enjoy from it with Python and PostGIS.

JTS has a minimum diameter method tsusiatsoftware.net/jts/javadoc/com/vividsolutions/jts/… that can also directly return "the minimum rectangular Polygon which encloses the input geometry". Convert you point set into multipoint and run the function. But based on your screen capture you may have a polygon to start with, and that you can use directly as input.

enter image description here

user30184
  • 65,331
  • 4
  • 65
  • 118
  • Didn't realize that. It's actually already in GEOS [1]. So, let's inform PostGIS developers about this...

    [1] http://geos.refractions.net/ro/doxygen_docs/html/classgeos_1_1algorithm_1_1MinimumDiameter.html

     – Stefan Feb 07 '15 at 14:53
1

I just was pointed to an answer by Rémi-C over at PostGIS mailinglist. He wrote:

"Here is a working solution in plpgsql, designed to be fast to code (and slow to execute :-/)"

https://github.com/Remi-C/PPPP_utilities/blob/6e9e8524812961b013b899466fe833dfa5d926e9/postgis/rc_oriented_bbox_deom_axis.sql

That function rc_BboxOrientedFromGeom(geom) does the job in PL/pgsql. Thanks to Rémi-C!

As he said and as I suggested above, most probably there exist more efficient solutions.

Stefan
  • 811
  • 7
  • 12
0

There's a cool solution for this: https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy

import numpy as np
from shapely.geometry.base import BaseGeometry
from shapely.geometry import Point, Polygon, MultiPolygon, GeometryCollection

def pca_eigenvectors(pts: np.ndarray) -> np.ndarray:
    """
    Returns the principal axes of a set of points.
    Method is essentially running a PCA on the points.


Parameters
----------
pts : array_like
"""
ca = np.cov(pts, y=None, rowvar=False, bias=True)
val, vect = np.linalg.eig(ca)

return np.transpose(vect)





def oriented_bounding_box(pts: np.ndarray) -> np.ndarray:
    """
    Returns the oriented bounding box width set of points.


Based on [Create the Oriented Bounding-box (OBB) with Python and NumPy](https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy).

Parameters
----------
pts : array_like
"""
tvect = pca_eigenvectors(pts)
rot_matrix = np.linalg.inv(tvect)

rot_arr = np.dot(pts, rot_matrix)

mina = np.min(rot_arr, axis=0)
maxa = np.max(rot_arr, axis=0)
diff = (maxa - mina) * 0.5

center = mina + diff

half_w, half_h = diff
corners = np.array([
    center + [-half_w, -half_h],
    center + [half_w, -half_h],
    center + [half_w, half_h],
    center + [-half_w, half_h],
])

return np.dot(corners, tvect)





def polygon_from_obb(obb: np.ndarray) -> Polygon:
    """
    Returns the oriented bounding box width set of points.


Parameters
----------
obb : array_like
"""
obb = np.vstack((obb, obb[0]))
return Polygon(obb)

This essentially performs a mini Principal Component Analysis to pull out the two perpendicular axes best describing the points, which can then be used to draw. It's all implemented with numpy for the most part so performant for most use cases.

I've pulled that example's implementation and added some utilities for working with shapely and geopandas here: https://github.com/raphaellaude/geo-obb/blob/main/geoobb/obb.py

See example usage: https://github.com/raphaellaude/geo-obb/blob/main/examples/parcel_obbs.ipynb

fishmulch
  • 116
  • 3