Question: 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.

Answer: 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:  sorted T based on whether P is above pj or not. 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 functions return 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()