用python擬合等角螺線的實現(xiàn)示例

r = fixed*np.exp(theta/np.tan(fixed))

import numpy as np
import matplotlib.pyplot as plt

from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['FangSong']
mpl.rcParams['axes.unicode_minus'] = False

def plotSpiral(core, fixed, phase=0, circle=4):
  """繪制等角螺線
  core		- 等角螺線的中心坐標(biāo)，tuple類型
  fixed    - 等角螺線的固定角度，單位：度（°）。fixed大于零則為順時針螺線，小于零則為逆時針螺線
  phase    - 初始相位，單位：圈（360°）。對順時針螺線，該數(shù)值越大，螺線越大，對逆時針螺線則相反
  circle   - 螺線可見部分的圈數(shù)，單位：圈（360°）
  """
  
  plt.axis("equal")
  plt.plot([core[0]], [core[1]], c='red', marker='+', markersize=10)
  
  fixed_rad = np.radians(90 + fixed)
  theta = np.linspace(0, circle*2*np.pi, 361) + phase*2*np.pi
  r = fixed_rad*np.exp(theta/np.tan(fixed_rad))
  x = r*np.cos(theta) + core[0]
  y = r*np.sin(theta) - core[1]
  plt.plot(x, y, c='blue')
  
  plt.show()

import numpy as np
from scipy import optimize

def fit_spiral(core, dots):
  """擬合等角螺線，返回定角fixed，初始相位phase"""
  
  fixed_ccw = 0.445*np.pi
  fixed_cw = 0.555*np.pi
  
  # 將dots拆分成x_list和y_list
  x_list, y_list = list(), list()
  for x, y in dots:
    x_list.append(x-core[0])
    y_list.append(y-core[1])
  
  # 計算距離
  x = np.array(x_list)
  y = np.array(y_list)
  r = np.hypot(x,y)
  
  # 按照距離排序
  sort_mask = np.argsort(r)
  x = x[sort_mask]
  y = y[sort_mask]
  r = r[sort_mask]
  
  # 計算角度
  theta = np.arctan(y/x)
  theta[x<0] += np.pi
  
  # 確定順序（CW-順時針，CCW-逆時針）
  d = np.diff(theta)
  print(d)
  ccw = d[d>0].size > d[d<0].size
  print('ccw=',ccw)
  
  # 調(diào)整角度為升序（CCW）或降序（CW）
  if ccw:
    for i in range(1, theta.size):
      while theta[i] < theta[i-1]:
        theta[i] += 2*np.pi
      
      dtheta = theta[i] - theta[i-1]
      while r[i]/r[i-1] > 1.8*np.exp(dtheta/np.tan(fixed_ccw)):
        theta[i] += 2*np.pi
        dtheta = theta[i] - theta[i-1]
  else:
    for i in range(theta.size-1)[::-1]:
      while theta[i] < theta[i+1]:
        theta[i] += 2*np.pi
      
      dtheta = theta[i+1] - theta[i]
      while r[i+1]/r[i] > 1.8*np.exp(dtheta/np.tan(fixed_cw)):
        theta[i] += 2*np.pi
        dtheta = theta[i+1] - theta[i]
  
  # 定義擬合函數(shù)
  def fmax(theta, fixed, phase):
    fixed = np.radians(90 + fixed)
    return fixed*np.exp((theta+phase*2*np.pi)/np.tan(fixed))
  
  try: 
    fita, fitb = optimize.curve_fit(fmax, theta, r, [2-int(ccw), 0], maxfev=10000)
    return fita
  except:
    return None

core = (530, 496)
dots = [(467,538), (448,675), (522,484), (513,451), (811,519)]
result = fit_spiral(core, dots)
if isinstance(result, np.ndarray):
  plotSpiral(core, result[0], phase=result[1], circle=4)
else:
  print(u'擬合失敗')