Question: I'm working with 2D shapes represented by their boundary contours (as ordered x,y coordinates), and I want to check if a shape is star-convex. If it is, I'd like to find a star center — i.e., a point from which the entire shape is visible (meaning: every line segment to every boundary point lies entirely within the shape).? Is there an elaborate non heuristic way of doing so? If there is already a python code for this then it would be even better :) This is what I have tried so far: I implemented a heuristic approach based on a midpoint visibility test. The idea is: A point is a valid star center if, for every boundary point, the midpoint of the line segment between them lies inside the polygon. I sample candidate points starting from the centroid, perturb them randomly within the bounding box, and check if one satisfies the condition. but this doesn't really guarantee finding the star-center or checking if it star-convex. There's also no guarantee that random sampling will find the right region (kernel) in the first place. def is_star_convex_and_get_center(x, y, n_samples=500): coords = np.column_stack((x, y)) poly = Polygon(coords) if not poly.is_valid or poly.area == 0: print("Invalid shape.") return None # Try centroid first candidates = [np.mean(coords, axis=0)] # Add some random points around the centroid bbox = poly.bounds scale = max(bbox[2] - bbox[0], bbox[3] - bbox[1]) for _ in range(n_samples - 1): offset = (np.random.rand(2) - 0.5) * scale candidates.append(candidates[0] + offset) # Check visibility condition for each candidate for pt in candidates: visible = True for px, py in coords: mid = 0.5 * (pt + np.array([px, py])) if not poly.contains(Point(mid)): visible = False break if visible: return pt # pt is a valid star center # No valid star center found return None 

Answer: I think you can do this by linear-programming.  Every CONVEX corner of your figure will give you a triangular wedge in which any solution must lie. In this wedge the X,Y point which is hopefully the centre must differ from the corner node by a combination of positive multiples of the two side vectors. i.e.  or  You then have to solve a matrix equation for X, Y and all the ak and bk, subject to all a's and b's being positive. This is exactly what scipy.optimize.linprog does: see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html In the following code IT IS ASSUMED THAT THE FIGURE IS TRACED ANTICLOCKWISE!!! import numpy as np from scipy.optimize import linprog import matplotlib.pyplot as plt def starConvex( x, y ): if x[0] != x[-1] or y[0] != y[-1]: # close the loop if necessary x.append( x[0] ) y.append( y[0] ) N = len( x ) - 1 # Count convex corners Nconvex = 0 for i in range( N ): i0, im, ip = i, i - 1, i + 1 if im < 0: im = N - 1 dxlow , dylow = x[im] - x[i0], y[im] - y[i0] dxhigh, dyhigh = x[ip] - x[i0], y[ip] - y[i0] if dxhigh * dylow - dxlow * dyhigh >= 0: Nconvex += 1 # Set up matrices for linear programming A = np.zeros( ( 2 * Nconvex, 2 * Nconvex + 2 ) ) b = np.zeros( 2 * Nconvex ) c = np.ones( 2 * Nconvex + 2 ); c[0] = c[1] = 0 k = 0 for i in range( N ): i0, im, ip = i, i - 1, i + 1 if im < 0: im = N - 1 dxlow , dylow = x[im] - x[i0], y[im] - y[i0] dxhigh, dyhigh = x[ip] - x[i0], y[ip] - y[i0] if dxhigh * dylow - dxlow * dyhigh >= 0: A[2*k ,0 ] = 1 A[2*k+1,1 ] = 1 A[2*k ,2*k+2] = -dxlow A[2*k ,2*k+3] = -dxhigh A[2*k+1,2*k+2] = -dylow A[2*k+1,2*k+3] = -dyhigh b[2*k ] = x[i0] b[2*k+1 ] = y[i0] k += 1 bounds=[(None,None),(None,None)] for _ in range( 2*Nconvex ): bounds.append( (0,None) ) res = linprog( c, A_eq=A, b_eq=b, bounds=bounds ) values = res.x if res.success: return True, values[0], values[1] else: return False, None, None xpts = [ 0.0, 1.0, 3.0, 1.0, 0.0, -1.0, -3.0, -1.0 ] ypts = [ -3.0, -1.0, 0.0, 1.0, 3.0, 1.0, 0.0, -1.0 ] # This succeeds #ypts = [ -3.0, -1.0, 0.0, 1.0, -1.0, 1.0, 0.0, -1.0 ] # This succeeds #ypts = [ -3.0, -1.0, 0.0, 1.0, -2.0, 1.0, 0.0, -1.0 ] # This fails test, x0, y0 = starConvex( xpts, ypts ) if test: print( "SUCCESS! x0, y0 = ", x0, y0 ) plt.plot( [x0], [y0], 'ro' ) else: print( "FAILURE" ) plt.plot( xpts, ypts, 'b' ) plt.show()  Here is the first case, which works OK (but you can't guarantee which feasible point will be chose as "star centre"; one possibility is shown in red.)  Below is the last case, where there is no point from which you can see the whole of the star.