python根据边界生成路径

如上图所示，根据边界生成路径的问题，不知道有没有大佬有好的数学方法。

我使用的方法是将两边的线向四边形内平移，再对得到的两边等距地取交点。这种方法对简单的图形可以很好地解决，复杂的图形很难处理（没试过）。使用的是python，代码在下面。刚学习python，求助帖，语法、格式有错误望见谅。

from scipy.linalg import solve
import math
import matplotlib.pyplot as plt
import numpy as np
putin=[[3,4], [5.5,3], [2,2], [4.5,1.2]]   #y值大的在前面
d=0.2

plt.plot([putin[0][0],putin[1][0]],[putin[0][1],putin[1][1]],c='r')
plt.plot([putin[0][0],putin[2][0]],[putin[0][1],putin[2][1]],c='r')
plt.plot([putin[3][0],putin[2][0]],[putin[3][1],putin[2][1]],c='r')
plt.plot([putin[3][0],putin[1][0]],[putin[3][1],putin[1][1]],c='r')



def down_down_pyhjiaodian(x0,y0,x1,y1,m0,n0,m1,n1,d):
    a11=(y1-y0)/(x1-x0)
    a12=-1
    b1=((y1-y0)/(x1-x0))*x0-y0+d*math.sqrt(((y1-y0)/(x1-x0))**2+1)
    a21=(n1-n0)/(m1-m0)
    a22=-1
    b2=((n1-n0)/(m1-m0))*m0-n0+d*math.sqrt(((n1-n0)/(m1-m0))**2+1)
    a=[[a11,a12],[a21,a22]]
    b=[b1,b2]
    x=solve(a,b)
    return x

def down_rise_pyhjiaodian(x0,y0,x1,y1,m0,n0,m1,n1,d):
    a11=(y1-y0)/(x1-x0)
    a12=-1
    b1=((y1-y0)/(x1-x0))*x0-y0+d*math.sqrt(((y1-y0)/(x1-x0))**2+1)
    a21=(n1-n0)/(m1-m0)
    a22=-1
    b2=((n1-n0)/(m1-m0))*m0-n0-d*math.sqrt(((n1-n0)/(m1-m0))**2+1)
    a=[[a11,a12],[a21,a22]]
    b=[b1,b2]
    x=solve(a,b)
    return x

def rise_down_pyhjiaodian(x0,y0,x1,y1,m0,n0,m1,n1,d):
    a11=(y1-y0)/(x1-x0)
    a12=-1
    b1=((y1-y0)/(x1-x0))*x0-y0-d*math.sqrt(((y1-y0)/(x1-x0))**2+1)
    a21=(n1-n0)/(m1-m0)
    a22=-1
    b2=((n1-n0)/(m1-m0))*m0-n0+d*math.sqrt(((n1-n0)/(m1-m0))**2+1)
    a=[[a11,a12],[a21,a22]]
    b=[b1,b2]
    x=solve(a,b)
    return x
def rise_rise_pyhjiaodian(x0,y0,x1,y1,m0,n0,m1,n1,d):
    a11=(y1-y0)/(x1-x0)
    a12=-1
    b1=((y1-y0)/(x1-x0))*x0-y0-d*math.sqrt(((y1-y0)/(x1-x0))**2+1)
    a21=(n1-n0)/(m1-m0)
    a22=-1
    b2=((n1-n0)/(m1-m0))*m0-n0-d*math.sqrt(((n1-n0)/(m1-m0))**2+1)
    a=[[a11,a12],[a21,a22]]
    b=[b1,b2]
    x=solve(a,b)
    return x



def point_line_distance(x0,y0,x1,y1,m,n):
    distance=(math.sqrt((((y1-y0)/(x1-x0))*m-n+y0-((y1-y0)/(x1-x0))*x0)**2))/(math.sqrt(((y1-y0)/(x1-x0))**2+1))
    return distance    

if (putin[0][1]-putin[2][1])/(putin[0][0]-putin[2][0])>0:
    dian3=down_rise_pyhjiaodian(putin[0][0],putin[0][1],putin[2][0],putin[2][1],putin[2][0],putin[2][1],putin[3][0],putin[3][1],d)
    dian1=down_down_pyhjiaodian(putin[0][0],putin[0][1],putin[2][0],putin[2][1],putin[0][0],putin[0][1],putin[1][0],putin[1][1],d)
else:
    dian3=rise_rise_pyhjiaodian(putin[0][0],putin[0][1],putin[2][0],putin[2][1],putin[2][0],putin[2][1],putin[3][0],putin[3][1],d)
    dian1=rise_down_pyhjiaodian(putin[0][0],putin[0][1],putin[2][0],putin[2][1],putin[0][0],putin[0][1],putin[1][0],putin[1][1],d)
if (putin[1][1]-putin[3][1])/(putin[1][0]-putin[3][0])<0:
    dian2=down_down_pyhjiaodian(putin[1][0],putin[1][1],putin[3][0],putin[3][1],putin[0][0],putin[0][1],putin[1][0],putin[1][1],d)
    dian4=down_rise_pyhjiaodian(putin[1][0],putin[1][1],putin[3][0],putin[3][1],putin[3][0],putin[3][1],putin[2][0],putin[2][1],d)
else:
    dian2=rise_down_pyhjiaodian(putin[1][0],putin[1][1],putin[3][0],putin[3][1],putin[0][0],putin[0][1],putin[1][0],putin[1][1],d)
    dian4=rise_rise_pyhjiaodian(putin[1][0],putin[1][1],putin[3][0],putin[3][1],putin[3][0],putin[3][1],putin[2][0],putin[2][1],d)


#print(dian1[0])
#print(dian2)
#print(dian3)
#print(dian4)
plt.scatter(dian1[0],dian1[1],c='r')
plt.scatter(dian2[0],dian2[1],c='r')
plt.scatter(dian3[0],dian3[1],c='r')
plt.scatter(dian4[0],dian4[1],c='r')

rel=[dian1,dian2]
l1=point_line_distance(dian3[0],dian3[1],dian4[0],dian4[1],dian1[0],dian1[1])
l2=point_line_distance(dian3[0],dian3[1],dian4[0],dian4[1],dian2[0],dian2[1])
i=1
#print(rel)
while l1>d and l2>d:
    #直线13
    f11=(dian1[1]-dian3[1])/(dian1[0]-dian3[0])
    f12=-1
    g1=((dian1[1]-dian3[1])/(dian1[0]-dian3[0]))*dian1[0]-dian1[1]
    #直线24
    f21=(dian2[1]-dian4[1])/(dian2[0]-dian4[0])
    f22=-1
    g2=((dian2[1]-dian4[1])/(dian2[0]-dian4[0]))*dian2[0]-dian2[1]
    #直线12
    f1=(dian1[1]-dian2[1])/(dian1[0]-dian2[0])
    f2=-1
    g=((dian1[1]-dian2[1])/(dian1[0]-dian2[0]))*dian1[0]-dian1[1]+2*d*math.sqrt(((dian1[1]-dian2[1])/(dian1[0]-dian2[0]))**2+1)
    
    A1=[[f11,f12],[f1,f2]]
    B1=[g1,g]
    
    A2=[[f21,f22],[f1,f2]]
    B2=[g2,g]
    
    dian1=solve(A1,B1)
    dian2=solve(A2,B2)
    
    
    if i%2==0:
        #print('偶数')
        rel.append(dian1)
        rel.append(dian2)
    else:
        rel.append(dian2)
        rel.append(dian1)
        
    #print('rel=',rel)
    l1=point_line_distance(dian3[0],dian3[1],dian4[0],dian4[1],dian1[0],dian1[1])
    l2=point_line_distance(dian3[0],dian3[1],dian4[0],dian4[1],dian2[0],dian2[1])
    i=i+1
    if dian1[1]<=dian3[1] or dian2[1]<=dian4[1]:
        break
    if dian1[1]<=dian4[1] or dian2[1]<=dian3[1]:
        break    

if i%2==0:
    #print('偶数')
    rel.append(dian3)
    rel.append(dian4) 
else:
    rel.append(dian4)
    rel.append(dian3)


rel= np.array(rel)  
print(rel)
print(rel[:,0])
print(rel[:,1])
print(i)
num=len(rel)
for i in range(num-1):
    plt.plot([rel[i][0],rel[i+1][0]],[rel[i][1],rel[i+1][1]])


plt.plot([dian3[0],dian4[0]],[dian3[1],dian4[1]])
plt.scatter(rel[:,0],rel[:,1],c='b')
plt.show()