Skip to content

Instantly share code, notes, and snippets.

@akabe

akabe/find2d_bl.py

Created January 28, 2025 05:23
Show Gist options
  • Save akabe/3cffc028e219262bfd2708a1fd92d832 to your computer and use it in GitHub Desktop.
Save akabe/3cffc028e219262bfd2708a1fd92d832 to your computer and use it in GitHub Desktop.
Download ZIP
「3 次元パッキングに対する効率的な bottom-left 法」の実装
Raw
find2d_bl.py
from __future__ import annotations
from dataclasses import dataclass
from enum import IntEnum
from math import ceil, log2
from typing import Optional
import pytest
@dataclass(frozen=True)
class Point2D:
x: int
y: int
@dataclass(frozen=True)
class Rect2D:
width: int
height: int
@dataclass(frozen=True)
class NFP2D:
x1: int
y1: int
x2: int
y2: int
@staticmethod
def create(point: Point2D, rect: Rect2D, new_rect: Rect2D) -> NFP2D:
x1 = point.x - new_rect.width
y1 = point.y - new_rect.height
x2 = point.x + rect.width
y2 = point.y + rect.height
return NFP2D(x1, y1, x2, y2)
class VEdgeType(IntEnum): # right < left になるように
RIGHT = 0
LEFT = 1
class HEdgeType(IntEnum): # top < bottom になるように
TOP = 0
BOTTOM = 1
@dataclass(frozen=True)
class VEdge:
x: int
y: int
type_: VEdgeType
placement_index: int
@dataclass(frozen=True)
class HEdge:
x: int
y: int
type_: HEdgeType
placement_index: int
## -----------------------------------------------------------------------------
##
## Find2D-BL で利用する完全二分木
##
## -----------------------------------------------------------------------------
@dataclass
class BLNode:
"""Find2D-BL で利用する完全二分木のノード"""
p_self: int
p_min: int
p_max: int
parent: Optional[BLNode]
left: Optional[BLNode]
right: Optional[BLNode]
def add(self, lambda_: int):
self.p_self += lambda_
self.p_min += lambda_
# self.p_max += lambda_ # 現論文には記載があるが、なくても BL 法は実装できる。
def update_min_max(self):
assert self.left is not None and self.right is not None
self.p_min = self.p_self + min(self.left.p_min, self.right.p_min)
# self.p_max = self.p_self + max(self.left.p_max, self.right.p_max)
class BLTree:
root: BLNode
leaves: list[BLNode]
def __init__(self, num_intervals: int):
"""Find2D-BL で利用する完全二分木の初期状態を生成する。
Parameters
----------
num_intervals
葉数
Returns
-------
初期状態の完全二分木
"""
depth = ceil(log2(num_intervals - 1))
num_leaves = 2**depth # 木が持つすべての葉の数
num_nodes = 2 ** (depth + 1) - 1 # 木が持つすべてのノード数
num_internal_nodes = num_nodes - num_leaves # 内部ノード数
num_dummy_leaves = num_leaves - num_intervals + 1 # ダミーノードの数
nodes = [BLNode(0, 0, 0, None, None, None) for _ in range(num_nodes)]
# 親子の参照を初期化する
for i in range(num_internal_nodes):
left = nodes[2 * i + 1]
right = nodes[2 * i + 2]
nodes[i].left = left
nodes[i].right = right
left.parent = nodes[i]
right.parent = nodes[i]
# ダミー葉は p_self=1 で初期化する（常に長方形を配置したくないノード）
for i in range(num_nodes - num_dummy_leaves, num_nodes):
nodes[i].p_self = 1
nodes[i].p_min = 1
nodes[i].p_max = 1
# 内部ノードは子の情報を元に値を決める
for i in reversed(range(num_internal_nodes)):
nodes[i].update_min_max()
self.root = nodes[0]
self.leaves = nodes[num_internal_nodes:]
def update_values(self, sweep_line_type: HEdgeType, left_index: int, right_index: int):
"""新たに遭遇した NFP の辺に応じて、葉ごとの NFP の重複数を更新する。
Parameters
----------
sweep_line_type
走査線のタイプ（上辺 or 下辺）
left_index
新たに遭遇した NFP の左辺に対応する葉のインデックス
right_index
新たに遭遇した NFP の右辺に対応する葉のインデックス
"""
assert left_index <= right_index
lambda_ = -1 if sweep_line_type == HEdgeType.TOP else +1
left = self.leaves[left_index]
right = self.leaves[right_index]
left.add(lambda_)
right.add(lambda_)
while left.parent is not right.parent:
assert left.parent is not None and right.parent is not None
assert left.parent.right is not None and right.parent.left is not None
# Add lambda_ to sibling nodes if necessary.
if left is not left.parent.right:
left.parent.right.add(lambda_)
if right is not right.parent.left:
right.parent.left.add(lambda_)
# Update parents.
left.parent.update_min_max()
right.parent.update_min_max()
left = left.parent
right = right.parent
# Update all nodes on the path to the root.
node = left.parent # Note that left.parent == right.parent
while node is not None:
node.update_min_max()
node = node.parent
def find_no_overwrap(self, alpha: int) -> Optional[int]:
"""左の葉から探索を始め、最初に NFP の重複数がゼロになる葉を返す。
Parameters
----------
alpha
探索範囲の左端の葉のインデックス
Returns
-------
alpha より右側で最初に NFP 重複数がゼロになる葉のインデックス。
そのような葉がない場合は、None を返す。"""
def _find(node: BLNode, offset: int, num_leaves: int) -> Optional[int]:
if alpha >= offset + num_leaves: # 探索範囲を外れた
return None
if node.p_min > 0: # node の配下に NFP の重複数が 0 の葉が存在しない
return None
if node.left is None: # node は葉
# p_min の定義より、葉について p_min == p_self なので node.p_self == 0 が保証される
return offset
else: # node は中間ノード
k = num_leaves // 2
# 左のノードの後で右のノードを探索する
result = _find(node.left, offset, k)
if result is None:
assert node.right is not None
result = _find(node.right, offset + k, k)
return result
return _find(self.root, 0, len(self.leaves))
## -----------------------------------------------------------------------------
##
## Find2D-BL のメインルーチン
##
## -----------------------------------------------------------------------------
def _NFPs_to_edges(nfps: list[NFP2D]) -> tuple[list[VEdge], list[HEdge]]:
"""NFP を辺に分解する。"""
left_right_edges = []
top_bottom_edges = []
for i, nfp in enumerate(nfps):
left = VEdge(nfp.x1, nfp.y2, VEdgeType.LEFT, i)
right = VEdge(nfp.x2, nfp.y2, VEdgeType.RIGHT, i)
top = HEdge(nfp.x2, nfp.y2, HEdgeType.TOP, i)
bottom = HEdge(nfp.x2, nfp.y1, HEdgeType.BOTTOM, i)
left_right_edges.extend([left, right])
top_bottom_edges.extend([top, bottom])
return left_right_edges, top_bottom_edges
def _get_interval_indices(intervals: list[VEdge]) -> tuple[list[int], list[int]]:
"""左辺・右辺の集合からインターバルを作る。"""
# i 番目の既配置図形 placements[i] に対して、
# - 左辺の x 座標は left_right_edges[ interval_indices[i][0] ]
# - 右辺の x 座標は left_right_edges[ interval_indices[i][1] ]
# になるように interval_indices を作る。
num_placements = len(intervals) // 2
left_indices = [-1 for _ in range(num_placements)]
right_indices = [-1 for _ in range(num_placements)]
for j, e in enumerate(intervals):
if e.type_ == VEdgeType.LEFT:
left_indices[e.placement_index] = j
else:
right_indices[e.placement_index] = j
return left_indices, right_indices
def find2d_bl(nfps: list[NFP2D], y_min: int = 0, y_max: int = 0) -> Optional[Point2D]:
"""Find2D-BL アルゴリズムを実行し、2 次元 BL 点を発見する。
Find2D-BL は既配置図形の NFP から BL 点を効率的に見つけるアルゴリズムである。
NFP の個数 n に対して、O(n log n) 時間で BL 点を発見できる。
Shinji Imahori, Yuyao Chien, Yuma Tanaka, and Mutsunori Yagiura.
Enumerating bottom-left stable positions for rectangle placements with overlap.
Journal of the Operations Research Society of Japan. Vol.57, No.1, March 2014, pp. 45--61.
https://www.jstage.jst.go.jp/article/jorsj/57/1/57_KJ00009350944/_pdf)
Parameters
----------
nfps
既配置図形と新規図形の NFP
y_max
探索範囲の y 座標の最大値
y_min
探索範囲の y 座標の最小値
"""
# NFP から辺を抽出
left_right_edges, top_bottom_edges = _NFPs_to_edges(nfps)
# 辺をソートする
intervals = sorted(left_right_edges, key=lambda e: (e.x, e.type_, e.y))
sweep_lines = sorted(top_bottom_edges, key=lambda e: (e.y, e.type_, e.x))
# sweep line の placement_index から左辺・右辺の index を引くための表を作る
placement_to_left_indices, placement_to_right_indices = _get_interval_indices(intervals)
# BLTree を作る
tree = BLTree(len(intervals))
# 下の辺から順に走査していく
for sweep_line in sweep_lines:
# 走査線に対応する既配置図形の index
i = sweep_line.placement_index
# 既配置図形の左辺・右辺に対応する葉ノードの index
left_leaf = placement_to_left_indices[i]
right_leaf = placement_to_right_indices[i] - 1
# NFP の重複数を更新する
tree.update_values(sweep_line.type_, left_leaf, right_leaf)
# 箱の中に新規図形を配置できなくなった
if sweep_line.y > y_max:
break
# 走査線がまだ箱の内側に入っていない
if sweep_line.y < y_min:
continue
# 既配置図形の上辺ならば、new_rect を配置可能かもしれない
if sweep_line.type_ == HEdgeType.TOP:
# NFP 重複がゼロになる点を探す
found_leaf = tree.find_no_overwrap(left_leaf)
if found_leaf is not None:
y = sweep_line.y
x = intervals[found_leaf].x
return Point2D(x, y)
return None # NFP 重複数がゼロになる点は見つからなかった
Raw
kawashima_11.py
from __future__ import annotations
import pytest
from bisect import insort
from typing import Iterable, Optional
from pydantic.dataclasses import dataclass
from find2d_bl import NFP2D, Point2D, Rect2D, find2d_bl
@dataclass(frozen=True)
class Point3D:
x: int
y: int
z: int
@dataclass(frozen=True)
class Rect3D:
width: int
depth: int
height: int
@dataclass(frozen=True)
class NFP3D:
x1: int
y1: int
z1: int
x2: int
y2: int
z2: int
@staticmethod
def create(point: Point3D, rect: Rect3D, new_rect: Rect3D) -> NFP3D:
x1 = point.x - new_rect.width
y1 = point.y - new_rect.depth
z1 = point.z - new_rect.height
x2 = point.x + rect.width
y2 = point.y + rect.depth
z2 = point.z + rect.height
return NFP3D(x1, y1, z1, x2, y2, z2)
def to_NFP2D(self) -> NFP2D:
return NFP2D(self.x1, self.y1, self.x2, self.y2)
@dataclass
class Face:
"""図形の上面・下面を表すクラス"""
x: int
y: int
z: int
placement_index: int
@staticmethod
def create_top_and_bottom(point: Point3D, rect: Rect3D, placement_index: int) -> tuple[Face, Face]:
z2 = point.z + rect.height
top_face = Face(
x=point.x,
y=point.y,
z=z2,
placement_index=placement_index,
)
bottom_face = Face(
x=point.x,
y=point.y,
z=point.z,
placement_index=placement_index,
)
return (top_face, bottom_face)
def bl_sort_key(self) -> tuple[int, int, int]:
return (self.z, self.y, self.x)
class DBLSolver:
box: Rect3D
placements: list[tuple[Point3D, Rect3D]]
top_faces: list[Face] # 既配置図形の上面（ソート済み）
bottom_faces: list[Face] # 既配置図形の下面（ソート済み）
def __init__(self, box: Rect3D):
self.box = box
self.placements = []
self.top_faces = []
self.bottom_faces = []
for point, rect in self.__create_placements_of_box(box):
self.put(point, rect)
@staticmethod
def __create_placements_of_box(box: Rect3D) -> list[tuple[Point3D, Rect3D]]:
L = 1 # 適当な厚み
w, d, h = box.width, box.depth, box.height
return [
(Point3D(-L, 0, 0), Rect3D(L, d, h)), # left
(Point3D(w, 0, 0), Rect3D(L, d, h)), # right
(Point3D(0, -L, 0), Rect3D(w, L, h)), # back
(Point3D(0, d, 0), Rect3D(w, L, h)), # front
(Point3D(0, 0, -L), Rect3D(w, d, L)), # bottom
(Point3D(0, 0, h), Rect3D(w, d, L)), # top
]
def put(self, point: Point3D, new_rect: Rect3D):
i = len(self.placements)
self.placements.append((point, new_rect))
top_face, bottom_face = Face.create_top_and_bottom(point, new_rect, i)
insort(
self.top_faces,
top_face,
key=Face.bl_sort_key,
)
insort(
self.bottom_faces,
bottom_face,
key=Face.bl_sort_key,
)
def find_BL_point(self, new_rect3d: Rect3D) -> Optional[Point3D]:
"""与えられた図形に対する BL 点を発見する。
Parameters
----------
new_rect3d
配置したい図形
Returns
-------
BL 点が存在するならば BL 点の座標を返す。存在しなければ None を返す。
"""
z_limit = self.box.depth - new_rect3d.depth
nfps = [NFP3D.create(point, rect, new_rect3d) for point, rect in self.placements]
E: set[NFP3D] = set()
it = 0 # the index of self.top_faces
ib = 0 # the index of self.bottom_faces
n = len(self.placements) # == len(self.top_faces) == len(self.bottom_faces)
while it < n and ib < n:
face_top = self.top_faces[it]
face_bot = self.bottom_faces[ib]
nfp_top = nfps[face_top.placement_index]
nfp_bot = nfps[face_bot.placement_index]
if nfp_bot.z1 < nfp_top.z2:
E.add(nfp_bot)
ib += 1
else:
E.remove(nfp_top)
# nfp_top の前面領域で Find2D-BL を実行
bl_point = self.__find2d_bl(new_rect3d, E, nfp_top, self.placements[face_top.placement_index])
# BL 点を発見できたら、結果を返す
if bl_point is not None:
return Point3D(bl_point.x, bl_point.y, face_top.z)
it += 1
return None
def __find2d_bl(
self, new_rect3d: Rect3D, nfps: Iterable[NFP3D], target_nfp: NFP3D, target_placement: tuple[Point3D, Rect3D]
) -> Optional[Point2D]:
# target_nfp と重なりのある NFP を抽出する。
overwrapped_nfps = self.__make_overwrapped_NFPs(nfps, target_nfp)
# target_nfp に対応する図形を箱とみなすので、箱の枠に対応する NFP を生成する。
new_rect2d = Rect2D(new_rect3d.width, new_rect3d.depth)
frames = self.__make_2d_frames(*target_placement)
frame_nfps = [NFP2D.create(point, rect, new_rect2d) for point, rect in frames]
# Find2D-BL の探索範囲を決める。
y_min = max(target_nfp.y1, 0)
y_max = min(target_nfp.y2, self.box.depth) - new_rect3d.depth
# Find2D-BL を実行する。
return find2d_bl(overwrapped_nfps + frame_nfps, y_min=y_min, y_max=y_max)
@staticmethod
def __make_2d_frames(point: Point3D, rect: Rect3D) -> list[tuple[Point2D, Rect2D]]:
x, y = point.x, point.y
w, d = rect.width, rect.depth
L = 1
return [
(Point2D(x - L, y), Rect2D(L, d)), # left frame
(Point2D(x + w, y), Rect2D(L, d)), # right frame
]
@staticmethod
def __make_overwrapped_NFPs(nfps: Iterable[NFP3D], target_nfp: NFP3D) -> list[NFP2D]:
overwrapped_nfps = []
for nfp in nfps:
x1 = max(nfp.x1, target_nfp.x1)
y1 = max(nfp.y1, target_nfp.y1)
x2 = min(nfp.x2, target_nfp.x2)
y2 = min(nfp.y2, target_nfp.y2)
if x1 < x2 and y1 < y2: # NFP の面積が 0 ではない
overwrapped_nfps.append(NFP2D(x1, y1, x2, y2))
return overwrapped_nfps
class TestDBLSolver:
@pytest.mark.parametrize(
["box", "placements", "new_rect", "expected"],
[
# 箱より十分小さい図形
pytest.param(
Rect3D(10, 10, 10),
[],
Rect3D(5, 5, 5),
Point3D(0, 0, 0),
id="rect_smaller_than_box",
),
# 箱と同じ大きさの図形
pytest.param(
Rect3D(10, 10, 10),
[],
Rect3D(10, 10, 10),
Point3D(0, 0, 0),
id="rect_of_the_same_size_as_box",
),
# 箱より width が大きい図形
pytest.param(
Rect3D(10, 10, 10),
[],
Rect3D(11, 5, 5),
None,
id="rect_with_width_larger_than_box",
),
# 箱より depth が大きい図形
pytest.param(
Rect3D(10, 10, 10),
[],
Rect3D(5, 11, 5),
None,
id="rect_with_depth_larger_than_box",
),
# 箱より height が大きい図形
pytest.param(
Rect3D(10, 10, 10),
[],
Rect3D(5, 5, 11),
None,
id="rect_with_height_larger_than_box",
),
# 既配置の図形の x 方向に配置
pytest.param(
Rect3D(10, 10, 10),
[
(Point3D(0, 0, 0), Rect3D(5, 5, 3)),
(Point3D(5, 0, 0), Rect3D(2, 3, 3)),
],
Rect3D(3, 4, 1),
Point3D(7, 0, 0),
id="place_in_x_direction_of_the_already_placed_rect",
),
# 既配置の図形の y 方向に配置
pytest.param(
Rect3D(10, 10, 10),
[
(Point3D(0, 0, 0), Rect3D(5, 5, 3)),
(Point3D(5, 0, 0), Rect3D(2, 3, 3)),
],
Rect3D(4, 4, 1),
Point3D(5, 3, 0),
id="place_in_y_direction_of_the_already_placed_rect",
),
# 既配置の図形の z 方向に配置
pytest.param(
Rect3D(10, 10, 10),
[
(Point3D(0, 0, 0), Rect3D(9, 5, 2)),
(Point3D(0, 5, 0), Rect3D(9, 5, 3)),
],
Rect3D(4, 4, 1),
Point3D(0, 0, 2),
id="place_in_z_direction_of_the_already_placed_rect",
),
# BL 点の x, y, z が全てゼロではない
pytest.param(
Rect3D(10, 10, 10),
[
(Point3D(0, 0, 0), Rect3D(2, 3, 3)),
(Point3D(2, 0, 0), Rect3D(6, 2, 4)),
(Point3D(2, 2, 0), Rect3D(5, 5, 1)),
],
Rect3D(4, 4, 2),
Point3D(2, 2, 1),
id="x_y_z_of_BL_point_are_not_zero",
),
# オーバーハングを許さない
pytest.param(
Rect3D(10, 10, 10),
[
(Point3D(0, 0, 0), Rect3D(5, 5, 3)),
],
Rect3D(9, 9, 1),
None,
id="prohibit_overhang",
),
# overwrapped_nfps の計算を間違えると失敗しやすいテストケース
pytest.param(
Rect3D(651, 411, 259),
[
(Point3D(0, 0, 0), Rect3D(330, 300, 240)),
(Point3D(330, 0, 0), Rect3D(300, 330, 40)),
(Point3D(330, 0, 40), Rect3D(300, 330, 20)),
(Point3D(0, 300, 0), Rect3D(180, 110, 100)),
(Point3D(180, 300, 0), Rect3D(90, 110, 100)),
],
Rect3D(90, 110, 100),
Point3D(x=330, y=0, z=60),
id="overwrapped_nfps_calculation",
),
],
)
def test_BL点を計算する(
self,
box: Rect3D,
placements: list[tuple[Point3D, Rect3D]],
new_rect: Rect3D,
expected: Optional[Point3D],
):
solver = DBLSolver(box)
for point, rect in placements:
solver.put(point, rect)
actual = solver.find_BL_point(new_rect)
assert expected == actual

Comments are disabled for this gist.