Compute minimum area convex k-gon in 2d

I am trying to solve the following problem: given a set of points P and a value k, find the area of the smallest convex k-gon defined by a subset of points S of P with |P| = n and |S| = k.

I found this paper describing an algorithm solving the problem in O(kn^3) (which is exactly what I need, considering in my case k can be big). Specifically, you can see the algorithm in section 4 (Algorithm 3).

Based on this paper, I have implemented the following python script:

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)

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))

# 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

I have also implemented a function computing the exact solution for P and k using a brute-force approach.

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

That being said, I have been stuck with this for a few days now. The solutions are sometimes correct and sometimes no. I guess the way I sort the points pl around pj is not correct (maybe?). I would really appreciate some help if any of you has experience with computational geometry or (even better) this paper/problem in particular.

I am providing a point set P and k for which the result of 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)
]

Thank you for your time and help!

P.S. you can expect P to be in general position and that only integer coordinates are allowed.

@Gianmarco Sorry for the late reply. I changed the function get_slope_clockwise(pj, P). The original function can not ensure the correct order using the absolute value of the slope. Absolute values can lead to points on opposite sides of the line (with positive and negative slopes) being treated as equivalent, which can disrupt the desired clockwise order. In the modified version:

1. sorted T based on whether P is above pj or not.
2. After that, sort based on the angle value as it is always in the range [0,pi].

def get_slope_clockwise(pj, P):
  T = [p for p in P if p != pj] #and p[0]>pj[0]
  #return sorted(T, key=lambda x: abs(slope(pj, x)))
  ref = np.array([-1, 0])
  #return sorted(T, key=lambda x: angle(ref, np.array(x) - np.array(pj)))
  #return sorted(T, key=lambda x: ((np.array(x) - np.array(pj))[1] < 0,angle(ref, np.array(x) - np.array(pj))))
  return sorted(
        T, 
        key=lambda x: (
            #(np.array(x) - np.array(pj))[1] < 0,  # Check if the point is below pj
            angle(ref, np.array(x)- np.array(pj)) if (np.array(x) - np.array(pj))[1] > 0 
            else 2*math.pi-angle( ref,np.array(x)- np.array(pj))  # Negate angle if below pj
        )
    )

Then,

import random
P = [(random.randrange(10**3), random.randrange(10**3)) for _ in P]
get_min_kgon(P,K) == reliable_min_kgon(P, K)[0]

In some cases, these two functions indeed returned different values, such as:

P=[(551, 728),
 (393, 966),
 (668, 453),
 (113, 616),
 (322, 33),
 (736, 665),
 (893, 440),
 (54, 597),
 (67, 140),
 (599, 614)]
K=5

Now, both function returned the same value in all local cases.

import matplotlib.pyplot as plt
result = reliable_min_kgon(P, K)
x, y = zip(*P)

# Plot the points
plt.figure(figsize=(8, 8))
plt.scatter(x, y, color='blue')

# Annotate the points
for i, (xi, yi) in enumerate(result[1]):
    plt.text(xi + 5, yi + 5, f'{i}', fontsize=12, color='red')

point0 = result[1][0]  # (551, 728)
point3 = result[1][3]  # (54, 597)
plt.plot([point0[0], point3[0]], [point0[1], point3[1]], color='green', linestyle='-', linewidth=2)
# Add grid, labels, and title
plt.grid(True)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scatter Plot of Points in P')

# Set equal scaling
plt.gca().set_aspect('equal', adjustable='box')

# Show the plot
plt.show()