Star-center of a star convex shape

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

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 a_k and b_k, 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.