Вычислить выпуклый k-угольник минимальной площади в 2d ⇐ Python

[Anonymous]

Я пытаюсь решить следующую задачу: учитывая набор точек P и значение k, найти площадь наименьшего выпуклого k-угольника, определенного подмножеством точек S из P с |P| = n и |S| = k.
Я нашел эту статью, описывающую алгоритм решения задачи за O(kn^3) (это именно то, что мне нужно, учитывая, что в моем случае k может быть большим). В частности, вы можете увидеть алгоритм в разделе 4 (Алгоритм 3).
На основе этой статьи я реализовал следующий скрипт Python:

Код: Выделить всё

import math

# Given a line defined by the points ab, it returns +1 if p is on the left side of the line, 0 if collinear and -1 if on the right side of the line
def side(a, b, p):
t = (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0])
return math.copysign(1, t)

# Computes the area of triangle abc
def area(a, b, c):
return abs(ConvexHull([a,b,c]).volume)

# Computes the slope of a line defined by points a,b
def slope(a, b):
return (b[1] - a[1]) / (b[0] - a[0])

# Returns the list of points in P above pi ordered clockwise around pi
def get_above_clockwise(pi, P, ref = np.array((-1, 0))):
above = [p for p in P if p[1] > pi[1]]
return sorted(above, key=lambda x: angle(ref, np.array(x) - np.array(pi)))

# Returns the list of points in P ordered clockwise (by slope) around pj
def get_slope_clockwise(pj, P):
T = [p for p in P if p != pj]
return sorted(T, key=lambda x: abs(slope(pj, x)))

def get_min_kgon(P, K):
D = {p:i for i, p in enumerate(P)}
T = np.full((K+1, len(P), len(P), len(P)), np.inf, dtype=np.float32)

total_min = np.inf

for pi in P:
T[2, D[pi]].flat = 0
PA = get_above_clockwise(pi, P)
for k in range(3, K+1):
for pj in PA:
min_area = np.inf

PA2 = get_slope_clockwise(pj, PA + [pi])
pi_idx = PA2.index(pi)
PA2 = PA2[(pi_idx+1):] + PA2[:pi_idx]

for pl in PA2:
if pl == pj: continue
if side(pi, pj, pl) == 1:
min_area = min(min_area, T[k-1, D[pi], D[pl], D[pj]] + area(pi, pj, pl))
T[k, D[pi], D[pj], D[pl]] = min_area

total_min = min(total_min, np.min(T[K, D[pi]].flat))

return total_min

Я также реализовал функцию, вычисляющую точное решение для P и k, используя метод грубой силы.

Код: Выделить всё

from itertools import combinations
from scipy.spatial import ConvexHull

def reliable_min_kgon(P, K):
m = np.inf
r = None
for lst in combinations(P, K):
ch = ConvexHull(lst)
if ch.volume < m and len(ch.vertices) == K:
m = ch.volume
r = lst
return m, r

При этом я застрял в этом уже несколько дней. Решения иногда правильные, иногда нет. Я думаю, что способ сортировки точек pl вокруг pj неверен (может быть?). Я был бы очень признателен за помощь, если кто-то из вас имеет опыт работы с вычислительной геометрией или (что еще лучше) с этой статьей/проблемой в частности.
Я предоставляю набор точек P и k, для которых результат get_min_kgon(P,k) != надежный_min_kgon(P, k).
P = [(102, 466), (435, 214), (860, 330), (270, 458), (106, 87), (71, 372), (700, 99), (20, 871), (614, 663), (121, 130)]
K = 5


Спасибо за ваше время и помощь!
P.S. вы можете ожидать, что P будет находиться в общем положении и разрешены только целочисленные координаты.

Подробнее здесь: https://stackoverflow.com/questions/788 ... -gon-in-2d