feature: Add target center, target ratio, and area tolerance to largest interior rectangle

largestinteriorrectangle/__init__.py

@@ -1,3 +1,3 @@
 from .lir import lir, pt1, pt2  # noqa: F401
 
-__version__ = "0.2.1"
+__version__ = "0.3.1"

largestinteriorrectangle/lir.py

@@ -3,7 +3,7 @@
 from .lir_within_polygon import largest_interior_rectangle as lir_within_polygon
 
 
-def lir(data, contour=None):
+def lir(data, contour=None, target_ratio=None, target_center=None, tolerance=None):
     """
     Computes the Largest Interior Rectangle.
     :param data: Can be
@@ -14,15 +14,39 @@ def lir(data, contour=None):
     :param contour: (optional) 2D ndarray with shape (n, 2) containing xy
     values of a specific contour where the rectangle could start (in all directions).
     Only needed for case 1.
+    :param target_ratio: (optional) float specifying the desired width/height ratio
+    of the rectangle. The rectangle with the largest area that has a width/height
+    ratio closest to the target_ratio is returned.
+    :param target_center: (optional) tuple of 2 floats specifying the desired center
+    of the rectangle. The rectangle with the largest area that has a center closest
+     to the target_center is returned.
+    :param tolerance: (optional) float specifying the tolerance for the target_center.
+    The tolerance with the largest area are considered.
     :return: 1D ndarray with lir specification: x, y, width, height
     :rtype: ndarray
     """
     if len(data.shape) == 3:
-        return lir_within_polygon(data)
+        return lir_within_polygon(
+            data,
+            target_ratio=target_ratio,
+            target_center=target_center,
+            tolerance=tolerance,
+        )
     if contour is None:
-        return lir_basis(data)
+        return lir_basis(
+            data,
+            target_ratio=target_ratio,
+            target_center=target_center,
+            tolerance=tolerance,
+        )
     else:
-        return lir_within_contour(data, contour)
+        return lir_within_contour(
+            data,
+            contour,
+            target_ratio=target_ratio,
+            target_center=target_center,
+            tolerance=tolerance,
+        )
 
 
 def pt1(lir):

largestinteriorrectangle/lir_basis.py

@@ -2,10 +2,20 @@
 import numpy as np
 
 
-def largest_interior_rectangle(grid):
+def largest_interior_rectangle(
+    grid, target_ratio=None, target_center=None, tolerance=None
+):
+    target_ratio = float(target_ratio) if target_ratio is not None else 0
     h_adjacency = horizontal_adjacency(grid)
     v_adjacency = vertical_adjacency(grid)
-    s_map = span_map(grid, h_adjacency, v_adjacency)
+    s_map = span_map(grid, h_adjacency, v_adjacency, target_ratio)
+    if target_center is not None:
+        assert (
+            tolerance is not None
+        ), "tolerance must be provided when target_center is used"
+        return biggest_span_in_span_map_closest_to_center(
+            s_map, target_center, tolerance
+        )
     return biggest_span_in_span_map(s_map)
 
 
@@ -77,6 +87,44 @@ def spans(h_vector, v_vector):
     return spans
 
 
+@nb.njit("uint32[:](uint32[:,:], float64)", cache=True)
+def biggest_constrained_span(spans, target_ratio):
+    if len(spans) == 0:
+        return np.array([0, 0], dtype=np.uint32)
+
+    max_area = 0
+    best_w = 0
+    best_h = 0
+
+    for i in nb.prange(len(spans)):
+        W_max = spans[i, 0]
+        H_max = spans[i, 1]
+
+        w1 = W_max
+        h1 = int(w1 / target_ratio)
+        if h1 > H_max:
+            h1 = H_max
+            w1 = int(h1 * target_ratio)
+
+        h2 = H_max
+        w2 = int(h2 * target_ratio)
+        if w2 > W_max:
+            w2 = W_max
+            h2 = int(w2 / target_ratio)
+
+        area1 = w1 * h1
+        area2 = w2 * h2
+
+        if area1 > max_area:
+            max_area = area1
+            best_w, best_h = w1, h1
+        if area2 > max_area:
+            max_area = area2
+            best_w, best_h = w2, h2
+
+    return np.array([best_w, best_h], dtype=np.uint32)
+
+
 @nb.njit("uint32[:](uint32[:,:])", cache=True)
 def biggest_span(spans):
     if len(spans) == 0:
@@ -87,11 +135,11 @@ def biggest_span(spans):
 
 
 @nb.njit(
-    "uint32[:, :, :](boolean[:,::1], uint32[:,::1], uint32[:,::1])",
+    "uint32[:, :, :](boolean[:,::1], uint32[:,::1], uint32[:,::1], float64)",
     parallel=True,
     cache=True,
 )
-def span_map(grid, h_adjacency, v_adjacency):
+def span_map(grid, h_adjacency, v_adjacency, target_ratio):
     y_values, x_values = grid.nonzero()
     span_map = np.zeros(grid.shape + (2,), dtype=np.uint32)
 
@@ -100,7 +148,10 @@ def span_map(grid, h_adjacency, v_adjacency):
         h_vec = h_vector(h_adjacency, x, y)
         v_vec = v_vector(v_adjacency, x, y)
         s = spans(h_vec, v_vec)
-        s = biggest_span(s)
+        if target_ratio > 0:
+            s = biggest_constrained_span(s, target_ratio)
+        else:
+            s = biggest_span(s)
         span_map[y, x, :] = s
 
     return span_map
@@ -114,3 +165,29 @@ def biggest_span_in_span_map(span_map):
     y = largest_rectangle_indices[0][0]
     span = span_map[y, x]
     return np.array([x, y, span[0], span[1]], dtype=np.uint32)
+
+
+def biggest_span_in_span_map_closest_to_center(span_map, target_center, tolerance=0.01):
+    areas = span_map[:, :, 0] * span_map[:, :, 1]
+
+    max_area = np.amax(areas)
+
+    if max_area == 0:
+        return np.array([0, 0, 0, 0], dtype=np.uint32)
+
+    threshold = max_area * (1.0 - tolerance)
+    ys, xs = np.where(areas >= threshold)
+
+    ws = span_map[ys, xs, 0]
+    hs = span_map[ys, xs, 1]
+
+    # inclusive rectangle convention (matches pt2)
+    rcx = xs.astype(np.float64) + (ws - 1.0) / 2.0
+    rcy = ys.astype(np.float64) + (hs - 1.0) / 2.0
+
+    cx, cy = float(target_center[0]), float(target_center[1])
+    d2 = (rcx - cx) ** 2 + (rcy - cy) ** 2
+
+    i = np.argmin(d2)
+
+    return np.array([xs[i], ys[i], ws[i], hs[i]], dtype=np.uint32)

largestinteriorrectangle/lir_within_contour.py

@@ -1,25 +1,45 @@
 import numba as nb
 import numpy as np
 
-from .lir_basis import biggest_span_in_span_map
+from .lir_basis import (
+    biggest_span_in_span_map,
+    biggest_span_in_span_map_closest_to_center,
+)
 from .lir_basis import h_vector as h_vector_top2bottom
 from .lir_basis import horizontal_adjacency as horizontal_adjacency_left2right
 from .lir_basis import predict_vector_size, span_map, spans
 from .lir_basis import v_vector as v_vector_left2right
 from .lir_basis import vertical_adjacency as vertical_adjacency_top2bottom
 
 
-def largest_interior_rectangle(grid, contour):
+def largest_interior_rectangle(
+    grid, contour, target_ratio=None, target_center=None, tolerance=None
+):
     adjacencies = adjacencies_all_directions(grid)
     contour = contour.astype("uint32", order="C")
 
-    s_map, _, saddle_candidates_map = create_maps(adjacencies, contour)
-    lir1 = biggest_span_in_span_map(s_map)
+    target_ratio = float(target_ratio) if target_ratio is not None else 0
+
+    s_map1, _, saddle_candidates_map = create_maps(adjacencies, contour, target_ratio)
+
+    s_map2 = span_map(
+        saddle_candidates_map, adjacencies[0], adjacencies[2], target_ratio
+    )
+
+    areas1 = s_map1[:, :, 0] * s_map1[:, :, 1]
+    areas2 = s_map2[:, :, 0] * s_map2[:, :, 1]
+    mask = areas2 > areas1
 
-    s_map = span_map(saddle_candidates_map, adjacencies[0], adjacencies[2])
-    lir2 = biggest_span_in_span_map(s_map)
+    s_map_combined = s_map1.copy()
+    s_map_combined[mask] = s_map2[mask]
+
+    if target_center is not None:
+        lir = biggest_span_in_span_map_closest_to_center(
+            s_map_combined, target_center, tolerance
+        )
+    else:
+        lir = biggest_span_in_span_map(s_map_combined)
 
-    lir = biggest_rectangle(lir1, lir2)
     return lir
 
 
@@ -123,28 +143,18 @@ def get_n_directions(spans_all_directions):
 
 
 @nb.njit(cache=True)
-def get_xy_array(x, y, spans, mode=0):
-    """0 - flip none, 1 - flip x, 2 - flip y, 3 - flip both"""
-    xy = spans.copy()
-    xy[:, 0] = x
-    xy[:, 1] = y
-    if mode == 1:
-        xy[:, 0] = xy[:, 0] - spans[:, 0] + 1
-    if mode == 2:
-        xy[:, 1] = xy[:, 1] - spans[:, 1] + 1
-    if mode == 3:
-        xy[:, 0] = xy[:, 0] - spans[:, 0] + 1
-        xy[:, 1] = xy[:, 1] - spans[:, 1] + 1
-    return xy
-
-
-@nb.njit(cache=True)
-def get_xy_arrays(x, y, spans_all_directions):
-    xy_l2r_t2b = get_xy_array(x, y, spans_all_directions[0], 0)
-    xy_r2l_t2b = get_xy_array(x, y, spans_all_directions[1], 1)
-    xy_l2r_b2t = get_xy_array(x, y, spans_all_directions[2], 2)
-    xy_r2l_b2t = get_xy_array(x, y, spans_all_directions[3], 3)
-    return xy_l2r_t2b, xy_r2l_t2b, xy_l2r_b2t, xy_r2l_b2t
+def get_top_left(x, y, w, h, direction_idx):
+    """0: none, 1: flip x, 2: flip y, 3: flip both"""
+    tx = x
+    ty = y
+    if direction_idx == 1:
+        tx = x - w + 1
+    elif direction_idx == 2:
+        ty = y - h + 1
+    elif direction_idx == 3:
+        tx = x - w + 1
+        ty = y - h + 1
+    return tx, ty
 
 
 @nb.njit(cache=True)
@@ -157,17 +167,22 @@ def cell_on_contour(x, y, contour):
 
 @nb.njit(
     "Tuple((uint32[:,:,::1], uint8[:,::1], boolean[:,::1]))"
-    "(UniTuple(uint32[:,::1], 4), uint32[:,::1])",
+    "(UniTuple(uint32[:,::1], 4), uint32[:,::1], float64)",
     parallel=True,
     cache=True,
 )
-def create_maps(adjacencies, contour):
+def create_maps(adjacencies, contour, target_ratio):
     h_left2right, h_right2left, v_top2bottom, v_bottom2top = adjacencies
 
     shape = h_left2right.shape
     span_map = np.zeros(shape + (2,), "uint32")
     direction_map = np.zeros(shape, "uint8")
     saddle_candidates_map = np.zeros(shape, "bool_")
+    constrained = target_ratio > 0
+
+    contour_grid = np.zeros(shape, "bool_")
+    for i in range(len(contour)):
+        contour_grid[contour[i, 1], contour[i, 0]] = True
 
     for idx in nb.prange(len(contour)):
         x, y = contour[idx, 0], contour[idx, 1]
@@ -176,17 +191,44 @@ def create_maps(adjacencies, contour):
         span_arrays = spans_all_directions(h_vectors, v_vectors)
         n = get_n_directions(span_arrays)
         direction_map[y, x] = n
-        xy_arrays = get_xy_arrays(x, y, span_arrays)
+
         for direction_idx in range(4):
-            xy_array = xy_arrays[direction_idx]
             span_array = span_arrays[direction_idx]
             for span_idx in range(span_array.shape[0]):
-                x, y = xy_array[span_idx][0], xy_array[span_idx][1]
-                w, h = span_array[span_idx][0], span_array[span_idx][1]
-                if w * h > span_map[y, x, 0] * span_map[y, x, 1]:
-                    span_map[y, x, :] = np.array([w, h], "uint32")
-                if n == 3 and not cell_on_contour(x, y, contour):
-                    saddle_candidates_map[y, x] = True
+                w_max, h_max = span_array[span_idx][0], span_array[span_idx][1]
+
+                if constrained:
+                    w1 = w_max
+                    h1 = int(w1 / target_ratio)
+                    if h1 > h_max:
+                        h1 = h_max
+                        w1 = int(h1 * target_ratio)
+
+                    h2 = h_max
+                    w2 = int(h2 * target_ratio)
+                    if w2 > w_max:
+                        w2 = w_max
+                        h2 = int(w2 / target_ratio)
+
+                    tx1, ty1 = get_top_left(x, y, w1, h1, direction_idx)
+                    if w1 * h1 > span_map[ty1, tx1, 0] * span_map[ty1, tx1, 1]:
+                        span_map[ty1, tx1, :] = np.array([w1, h1], "uint32")
+
+                    tx2, ty2 = get_top_left(x, y, w2, h2, direction_idx)
+                    if w2 * h2 > span_map[ty2, tx2, 0] * span_map[ty2, tx2, 1]:
+                        span_map[ty2, tx2, :] = np.array([w2, h2], "uint32")
+                else:
+                    w, h = w_max, h_max
+                    tx, ty = get_top_left(x, y, w, h, direction_idx)
+                    if w * h > span_map[ty, tx, 0] * span_map[ty, tx, 1]:
+                        span_map[ty, tx, :] = np.array([w, h], "uint32")
+
+            for dy in range(-1, 2):
+                for dx in range(-1, 2):
+                    ny, nx = y + dy, x + dx
+                    if 0 <= ny < shape[0] and 0 <= nx < shape[1]:
+                        if h_left2right[ny, nx] > 0 and not contour_grid[ny, nx]:
+                            saddle_candidates_map[ny, nx] = True
 
     return span_map, direction_map, saddle_candidates_map
 
@@ -197,3 +239,29 @@ def biggest_rectangle(*args):
         if rect[2] * rect[3] > biggest_rect[2] * biggest_rect[3]:
             biggest_rect = rect
     return biggest_rect
+
+
+def rectangle_center_distance(rect, target_center):
+    rcx = float(rect[0]) + (float(rect[2]) - 1.0) / 2.0
+    rcy = float(rect[1]) + (float(rect[3]) - 1.0) / 2.0
+    dx = rcx - float(target_center[0])
+    dy = rcy - float(target_center[1])
+    return dx * dx + dy * dy
+
+
+def choose_rectangle(rect1, rect2, target_center=None):
+    if target_center is None:
+        return biggest_rectangle(rect1, rect2)
+
+    d1 = rectangle_center_distance(rect1, target_center)
+    d2 = rectangle_center_distance(rect2, target_center)
+    if d1 < d2:
+        return rect1
+    if d2 < d1:
+        return rect2
+
+    a1 = rect1[2] * rect1[3]
+    a2 = rect2[2] * rect2[3]
+    if a1 >= a2:
+        return rect1
+    return rect2