二维最小面积凸k边形算法实现异常求助

我正在解决的问题：给定点集P与数值k，求由P的k点子集S构成的最小面积凸k边形的面积（其中|P|=n，|S|=k）。

我找到一篇论文，其描述的算法可在O(kn³)复杂度下解决该问题（符合我k值较大的场景需求），具体算法见第4节（Algorithm 3）。基于该论文，我实现了如下Python脚本：

import math
import numpy as np
from scipy.spatial import ConvexHull

# 判断点p在直线ab的哪一侧：+1为左侧，0为共线，-1为右侧
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)

def angle(v1, v2):
    v1_u = v1 / np.linalg.norm(v1)
    v2_u = v2 / np.linalg.norm(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

# 计算三角形abc的面积
def area(a, b, c):
    return abs(ConvexHull([a,b,c]).volume)

# 计算直线ab的斜率
def slope(a, b):
    return (b[1] - a[1]) / (b[0] - a[0])

# 返回P中位于pi上方、围绕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)))

# 返回围绕pj按斜率顺时针排序的P中点列表（排除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

我还实现了暴力求解函数以验证结果：

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

目前我已卡壳数日，代码结果时而正确时而错误，我猜测是围绕pj排序点pl的逻辑存在问题，恳请有计算几何经验或了解该论文/问题的人士提供帮助。

提供一组测试用例，其中get_min_kgon(P,k) != reliable_min_kgon(P, k)：

K = 5
P = [
(102, 466), (435, 214), (860, 330), (270, 458), (106, 87),
(71, 372), (700, 99), (20, 871), (614, 663), (121, 130)
]

注：输入点集处于一般位置，且仅包含整数坐标。

内容的提问来源于stack exchange，提问作者Toby