[Tex/LaTex] How to generate random closed shape and draw a line bisecting it

pgfplotstikz-pgf

Similar to this question, except I also want a single line bisecting the shape in both area and perimeter. So the dream is to end up with a shape like the ones shown here, and a single line through the shape, where the area of the shape to one side of the line is equal to the other and ditto for perimeter.

I don't even know if this is possible or not. How would one even begin to compute the area/perimeter of the shape?

Best Answer

I think this is a case where it might be wise to use an external tool for calculating the coordinates.

You can use Pythons shapely library together with the scipy.optimize tools for finding a split line:

import numpy as np
from shapely.geometry import Polygon, LineString
from shapely.ops import polygonize
import scipy.optimize

### Gaussian smoother from http://www.swharden.com/blog/2008-11-17-linear-data-smoothing-in-python/ for getting a nice polygon

def smoothListGaussian(list,degree=5):  
     window=degree*2-1  
     weight=np.array([1.0]*window)  
     weightGauss=[]  
     for i in range(window):  
         i=i-degree+1  
         frac=i/float(window)  
         gauss=1/(np.exp((4*(frac))**2))  
         weightGauss.append(gauss)  
     weight=np.array(weightGauss)*weight  
     smoothed=[0.0]*(len(list)-window)  
     for i in range(len(smoothed)):  
         smoothed[i]=sum(np.array(list[i:i+window])*weight)/sum(weight)  
     return smoothed  

# Generate the polygon
theta = np.linspace(0,2*np.pi,200, endpoint=False)
r = np.random.lognormal(0,0.4,200)
r = np.pad(r,(9,10),mode='wrap')

r = smoothListGaussian(r, degree=10)

coords = zip(np.cos(theta)*r, np.sin(theta)*r)

polygon = Polygon(coords)


# The function for splitting the polygon and calculating the objective function value
def splitPoly(poly, parameters):
    rot = parameters[0]
    shift = parameters[1]

    c = poly.centroid.coords[0]
    l = polygon.length
    shiftvec = (np.cos(rot), np.sin(rot))

    linestart = (c[0] - shiftvec[0]*l + shiftvec[1]*shift, c[1] - shiftvec[1]*l - shiftvec[0]*shift)
    lineend = (c[0] + shiftvec[0]*l + shiftvec[1]*shift, c[1] + shiftvec[1]*l - shiftvec[0]*shift)

    line = LineString([linestart, lineend])

    splitpoly = list(polygonize(poly.boundary.union( line ) ))
    areadiff = 1+np.abs(splitpoly[0].area - splitpoly[1].area)
    lengthdiff = 1+np.abs(splitpoly[0].length - splitpoly[1].length)
    return areadiff*lengthdiff-1, splitpoly

def wrapper(parameters):
    return splitPoly(polygon, parameters)[0]

# Use a grid search for finding a good starting value
res = scipy.optimize.brute(wrapper, ((0.0,2),(-0.1,0.1)))

# Nelder-Mead for optimizing the parameters
res = scipy.optimize.minimize(wrapper, res, method='nelder-mead')

# Split the polygon using the final parameter values
value, splitpoly = splitPoly(polygon, res.x)

print "Areas: ", splitpoly[0].area, splitpoly[1].area
print "Perimeters: ", splitpoly[0].length, splitpoly[1].length

# Write the coordinates
np.savetxt('polygonA.txt', np.array(list(splitpoly[0].exterior.coords)))
np.savetxt('polygonB.txt', np.array(list(splitpoly[1].exterior.coords)))

The polygons can then be plotted using PGFPlots:

\documentclass{standalone}
\usepackage{pgfplots}

\begin{document}
\begin{tikzpicture}
\begin{axis}[axis equal image, hide axis]
\addplot [no markers, fill=yellow] table {polygonA.txt};
\addplot [no markers, fill=orange] table {polygonB.txt};
\end{axis}
\end{tikzpicture}
\end{document}

Related Question