diff --git a/largestinteriorrectangle/lir_basis.py b/largestinteriorrectangle/lir_basis.py index 6eea98f..5e0239d 100644 --- a/largestinteriorrectangle/lir_basis.py +++ b/largestinteriorrectangle/lir_basis.py @@ -1,116 +1,169 @@ -import numba as nb +import math + import numpy as np -def largest_interior_rectangle(grid): - h_adjacency = horizontal_adjacency(grid) - v_adjacency = vertical_adjacency(grid) - s_map = span_map(grid, h_adjacency, v_adjacency) - return biggest_span_in_span_map(s_map) - - -@nb.njit("uint32[:,::1](boolean[:,::1])", parallel=True, cache=True) -def horizontal_adjacency(grid): - result = np.zeros((grid.shape[0], grid.shape[1]), dtype=np.uint32) - for y in nb.prange(grid.shape[0]): - span = 0 - for x in range(grid.shape[1] - 1, -1, -1): - if grid[y, x]: - span += 1 - else: - span = 0 - result[y, x] = span - return result - - -@nb.njit("uint32[:,::1](boolean[:,::1])", parallel=True, cache=True) -def vertical_adjacency(grid): - result = np.zeros((grid.shape[0], grid.shape[1]), dtype=np.uint32) - for x in nb.prange(grid.shape[1]): - span = 0 - for y in range(grid.shape[0] - 1, -1, -1): - if grid[y, x]: - span += 1 - else: - span = 0 - result[y, x] = span - return result - - -@nb.njit("uint32(uint32[:])", cache=True) -def predict_vector_size(array): - zero_indices = np.where(array == 0)[0] - if len(zero_indices) == 0: - if len(array) == 0: - return 0 - return len(array) - return zero_indices[0] - - -@nb.njit("uint32[:](uint32[:,::1], uint32, uint32)", cache=True) -def h_vector(h_adjacency, x, y): - vector_size = predict_vector_size(h_adjacency[y:, x]) - h_vector = np.zeros(vector_size, dtype=np.uint32) - h = np.inf - for p in range(vector_size): - h = np.minimum(h_adjacency[y + p, x], h) - h_vector[p] = h - h_vector = np.unique(h_vector)[::-1] - return h_vector - - -@nb.njit("uint32[:](uint32[:,::1], uint32, uint32)", cache=True) -def v_vector(v_adjacency, x, y): - vector_size = predict_vector_size(v_adjacency[y, x:]) - v_vector = np.zeros(vector_size, dtype=np.uint32) - v = np.inf - for q in range(vector_size): - v = np.minimum(v_adjacency[y, x + q], v) - v_vector[q] = v - v_vector = np.unique(v_vector)[::-1] - return v_vector - - -@nb.njit("uint32[:,:](uint32[:], uint32[:])", cache=True) -def spans(h_vector, v_vector): - spans = np.stack((h_vector, v_vector[::-1]), axis=1) - return spans - - -@nb.njit("uint32[:](uint32[:,:])", cache=True) -def biggest_span(spans): - if len(spans) == 0: - return np.array([0, 0], dtype=np.uint32) - areas = spans[:, 0] * spans[:, 1] - biggest_span_index = np.where(areas == np.amax(areas))[0][0] - return spans[biggest_span_index] - - -@nb.njit( - "uint32[:, :, :](boolean[:,::1], uint32[:,::1], uint32[:,::1])", - parallel=True, - cache=True, -) -def span_map(grid, h_adjacency, v_adjacency): - y_values, x_values = grid.nonzero() - span_map = np.zeros(grid.shape + (2,), dtype=np.uint32) - - for idx in nb.prange(len(x_values)): - x, y = x_values[idx], y_values[idx] - 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) - span_map[y, x, :] = s - - return span_map - - -@nb.njit("uint32[:](uint32[:, :, :])", cache=True) -def biggest_span_in_span_map(span_map): - areas = span_map[:, :, 0] * span_map[:, :, 1] - largest_rectangle_indices = np.where(areas == np.amax(areas)) - x = largest_rectangle_indices[1][0] - y = largest_rectangle_indices[0][0] - span = span_map[y, x] - return np.array([x, y, span[0], span[1]], dtype=np.uint32) +def _trim_grid(grid: np.array, center, target_ratio=0.00): + assert target_ratio >= 0, "target_ratio cannot be negative!" + rows, cols = grid.shape + center_x, center_y = center + + left_edge = center_x + for x in range(center_x, -1, -1): + if not grid[center_y, x]: + break + left_edge = x + + right_edge = center_x + for x in range(center_x, cols): + if not grid[center_y, x]: + break + right_edge = x + + up_edge = center_y + for y in range(center_y, -1, -1): + if not grid[y, center_x]: + break + up_edge = y + + down_edge = center_y + for y in range(center_y, rows): + if not grid[y, center_x]: + break + down_edge = y + + half_trimmed_x_range = min(right_edge - center_x, center_x - left_edge) + half_trimmed_y_range = min(down_edge - center_y, center_y - up_edge) + x_min = max(center_x - half_trimmed_x_range, 0) + x_max = min(center_x + half_trimmed_x_range, cols - 1) + y_min = max(center_y - half_trimmed_y_range, 0) + y_max = min(center_y + half_trimmed_y_range, rows - 1) + + if target_ratio: + w = x_max - x_min + h = y_max - y_min + if w > h: + w = h * target_ratio + x_min = int(center_x - w // 2) + x_max = int(center_x + w // 2) + else: + h = w / target_ratio + y_min = int(center_y - h // 2) + y_max = int(center_y + h // 2) + + trimmed_grid = np.zeros_like(grid) + trimmed_grid[y_min : y_max + 1, x_min : x_max + 1] = grid[ + y_min : y_max + 1, x_min : x_max + 1 + ] + return trimmed_grid + + +def _get_step(target_ratio, step_max=20): + scale_factor = 100000 # Scale up to avoid floating point precision issues + x_step = y_step = step_max + 1 + while (x_step > step_max or y_step > step_max) and scale_factor: + width = target_ratio + height = int(scale_factor) + width = int(width * scale_factor) + width = width + 1 if width % 2 else width + + gcd = math.gcd(width, height) # Compute the GCD + + x_step = width // gcd + y_step = height // gcd + scale_factor /= 10 + + return x_step, y_step + + +def largest_interior_rectangle( + grid: np.array, target_ratio=0.00, ratio_tolerance=0.01, target_center_coord=None +) -> np.array: + rows, cols = grid.shape + max_area = 0 + max_rect = (0, 0, 0, 0) # (x1, y1, x2, y2) + + if target_center_coord: + # Check if center_coord is in the grid + center_x, center_y = target_center_coord + if ( + 0 <= center_x < cols + and 0 <= center_y < rows + and grid[center_y, center_x] == 1 + ): + grid = _trim_grid(grid, target_center_coord, target_ratio=target_ratio) + x_step, y_step = _get_step(target_ratio) + x_extend = y_extend = 0 + # TODO: simplify this ugly while and that ugly if in it + while ( + center_y + y_extend < rows + and center_y - y_extend >= 0 + and center_x + x_extend < cols + and center_x - x_extend >= 0 + and grid[center_y + y_extend, center_x] + and grid[center_y - y_extend, center_x] + and grid[center_y, center_x + x_extend] + and grid[center_y, center_x - x_extend] + ): + y_extend += y_step + x_extend += x_step + if ( + center_y - y_extend >= 0 + and center_y + y_extend < rows + and center_x - x_extend >= 0 + and center_x + x_extend < cols + and grid[ + center_y - y_extend : center_y + y_extend + 1, + center_x - x_extend : center_x + x_extend + 1, + ].all() + ): + max_area = max(max_area, (x_extend * 2) * (y_extend * 2)) + max_rect = ( + center_x - x_extend, + center_y - y_extend, + center_x + x_extend, + center_y + y_extend, + ) + else: + print(f"⚠️Center coordinate is not in the grid! ({grid.shape = })") + + else: # auto mode + # ref: LeetCode 84. Largest Rectangle in Histogram + heights = np.zeros(cols) + for y2 in range(rows): + if ( + rows - y2 + ) * cols < max_area: # The remaining area is smaller than the max area + break + + heights = (heights + 1) * grid[y2, :] # compute the heights (histogram) + stack = [] # reset stack at every row + for x2 in range(cols + 1): + if (rows - y2) * ( + cols + 1 - x2 + ) < max_area: # The remaining area is smaller than the max area + break + + while stack and (x2 == cols or heights[stack[-1]] > heights[x2]): + height = heights[stack.pop()] + width = x2 if not stack else x2 - stack[-1] - 1 + curr_ratio = ( + width / height if height > 0 else 0 + ) # prevent ZeroDivisionError + if ( + target_ratio * (1 - ratio_tolerance) + <= curr_ratio + <= target_ratio * (1 + ratio_tolerance) + or target_ratio == 0 + ): + area = width * height + x1 = 0 if not stack else stack[-1] + 1 + rect = (x1, y2 - height + 1, x2 - 1, y2) + if area > max_area: + max_area = area + max_rect = rect + + stack.append(x2) + + x1, y1, x2, y2 = max_rect + return np.array([x1, y1, x2 - x1, y2 - y1], dtype="int")