Аппроксимация данных с помощью многосегментной кубической кривой Безье и расстояния, а также ограничения кривизны

Я нашел решение, которое соответствует моим критериям. Решение состоит в том, чтобы сначала найти B-сплайн, который аппроксимирует точки по методу наименьших квадратов, а затем преобразовать этот сплайн в многосегментную кривую Безье. B-сплайны имеют то преимущество, что в отличие от кривых Безье они не проходят через контрольные точки, а также предоставляют способ задать желаемую «гладкость» аппроксимационной кривой. Необходимая функциональность для создания такого сплайна реализована в библиотеке FITPACK, к которой scipy предлагает привязку python. Предположим, я прочитал свои данные в списках x и y, тогда я могу:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
tck,u = interpolate.splprep([x,y],s=3)
unew = np.arange(0,1.01,0.01)
out = interpolate.splev(unew,tck)
plt.figure()
plt.plot(x,y,out[0],out[1])
plt.show()

Результат выглядит так:

Если я хочу, чтобы кривая была более плавной, я могу увеличить параметр s до splprep. Если я хочу, чтобы приближение было ближе к данным, я могу уменьшить параметр s для меньшей плавности. Программно перебирая несколько s параметров, я могу найти хороший параметр, который соответствует заданным требованиям.

Однако вопрос в том, как преобразовать этот результат в кривую Безье. Ответ в этом письме Захари Пинкуса. Я воспроизведу его решение здесь, чтобы дать исчерпывающий ответ на мой вопрос:

def b_spline_to_bezier_series(tck, per = False):
  """Convert a parametric b-spline into a sequence of Bezier curves of the same degree.

  Inputs:
    tck : (t,c,k) tuple of b-spline knots, coefficients, and degree returned by splprep.
    per : if tck was created as a periodic spline, per *must* be true, else per *must* be false.

  Output:
    A list of Bezier curves of degree k that is equivalent to the input spline. 
    Each Bezier curve is an array of shape (k+1,d) where d is the dimension of the
    space; thus the curve includes the starting point, the k-1 internal control 
    points, and the endpoint, where each point is of d dimensions.
  """
  from fitpack import insert
  from numpy import asarray, unique, split, sum
  t,c,k = tck
  t = asarray(t)
  try:
    c[0][0]
  except:
    # I can't figure out a simple way to convert nonparametric splines to 
    # parametric splines. Oh well.
    raise TypeError("Only parametric b-splines are supported.")
  new_tck = tck
  if per:
    # ignore the leading and trailing k knots that exist to enforce periodicity 
    knots_to_consider = unique(t[k:-k])
  else:
    # the first and last k+1 knots are identical in the non-periodic case, so
    # no need to consider them when increasing the knot multiplicities below
    knots_to_consider = unique(t[k+1:-k-1])
  # For each unique knot, bring it's multiplicity up to the next multiple of k+1
  # This removes all continuity constraints between each of the original knots, 
  # creating a set of independent Bezier curves.
  desired_multiplicity = k+1
  for x in knots_to_consider:
    current_multiplicity = sum(t == x)
    remainder = current_multiplicity%desired_multiplicity
    if remainder != 0:
      # add enough knots to bring the current multiplicity up to the desired multiplicity
      number_to_insert = desired_multiplicity - remainder
      new_tck = insert(x, new_tck, number_to_insert, per)
  tt,cc,kk = new_tck
  # strip off the last k+1 knots, as they are redundant after knot insertion
  bezier_points = numpy.transpose(cc)[:-desired_multiplicity]
  if per:
    # again, ignore the leading and trailing k knots
    bezier_points = bezier_points[k:-k]
  # group the points into the desired bezier curves
  return split(bezier_points, len(bezier_points) / desired_multiplicity, axis = 0)

Так что B-Splines, FITPACK, numpy и scipy спасли мне день :)

1. полигонизировать данные

   найдите порядок точек, чтобы вы просто находили самые близкие друг к другу точки и пытались соединить их «линиями». Избегайте возврата к исходной точке

2. вычислить вывод по пути

   это изменение направления «линий», где вы нажимаете локальный минимум или максимум, это ваша контрольная точка ... Сделайте это, чтобы уменьшить ваши входные данные (оставьте только контрольные точки).

3. кривая

   теперь используйте эти точки как контрольные. Я настоятельно рекомендую использовать полином интерполяции как для x, так и для y по отдельности, например, примерно так:

   x=a0+a1*t+a2*t*t+a3*t*t*t
   y=b0+b1*t+b2*t*t+b3*t*t*t
   

   где a0..a3 вычисляются следующим образом:

   d1=0.5*(p2.x-p0.x);
   d2=0.5*(p3.x-p1.x);
   a0=p1.x;
   a1=d1;
   a2=(3.0*(p2.x-p1.x))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-p2.x+p1.x));
   

   • b0 .. b3 are computed in same way but use y coordinates of course
   • p0..p3 - контрольные точки для кривой кубической интерполяции
   • t =<0.0,1.0> - параметр кривой от p1 до p2

   это гарантирует, что позиция и первая деривация будут непрерывными (c1), а также вы можете использовать BEZIER, но это будет не так хорошо, как это.

[edit1] слишком острые края - БОЛЬШАЯ проблема

Чтобы решить эту проблему, вы можете удалить точки из набора данных до получения контрольных точек. Я могу придумать два способа сделать это прямо сейчас ... выберите, что лучше для вас

1. удалить точки из набора данных со слишком высокой первой производной

   dx/dl или dy/dl, где x,y - координаты, а l - длина кривой (вдоль ее пути). Точное вычисление радиуса кривизны на основе построения кривой сложно.

2. удалить точки из набора данных, которые приводят к слишком маленькому радиусу кривизны

   вычислить пересечение средней точки соседних сегментов линии (черные линии). Перпендикулярные оси, как на изображении (красные линии), расстояние между ними и точкой соединения (синяя линия) - это ваш радиус кривизны. Когда радиус кривизны меньше, то ваш предел удалите эту точку ...

   

   теперь, если вам действительно нужны только кубики BEZIER, вы можете преобразовать мою кубическую интерполяцию в кубическую BEZIER следующим образом:

//  ---------------------------------------------------------------------------
//  x=cx[0]+(t*cx[1])+(tt*cx[2])+(ttt*cx[3]); // cubic x=f(t), t = <0,1>
//  ---------------------------------------------------------------------------
//  cubic matrix                           bz4 = it4
//  ---------------------------------------------------------------------------
//  cx[0]=                            (    x0) =                    (    X1)
//  cx[1]=                   (3.0*x1)-(3.0*x0) =           (0.5*X2)         -(0.5*X0)
//  cx[2]=          (3.0*x2)-(6.0*x1)+(3.0*x0) = -(0.5*X3)+(2.0*X2)-(2.5*X1)+(    X0)
//  cx[3]= (    x3)-(3.0*x2)+(3.0*x1)-(    x0) =  (0.5*X3)-(1.5*X2)+(1.5*X1)-(0.5*X0)
//  ---------------------------------------------------------------------------
    const double m=1.0/6.0;
    double x0,y0,x1,y1,x2,y2,x3,y3;
    x0 = X1;           y0 = Y1;
    x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
    x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
    x3 = X2;           y3 = Y2;

Если вам нужно обратное преобразование, см .:

• Кривая Безье с контрольными точками внутри кривой

Вопрос был опубликован давно, но вот простое решение, основанное на splprep, нахождение минимального значения s, позволяющего выполнить критерии минимального радиуса кривизны.

route - это набор входных точек, первое измерение - это количество точек.

import numpy as np
from scipy.interpolate import splprep, splev

#The minimum curvature radius we want to enforce
minCurvatureConstraint = 2000
#Relative tolerance on the radius
relTol = 1.e-6

#Initial values for bisection search, should bound the solution
s_0 = 0
minCurvature_0 = 0
s_1 = 100000000 #Should be high enough to produce curvature radius larger than constraint
s_1 *= 2
minCurvature_1 = np.float('inf')

while np.abs(minCurvature_0 - minCurvature_1)>minCurvatureConstraint*relTol:
    s = 0.5 * (s_0 + s_1)
    tck, u  = splprep(np.transpose(route), s=s)
    smoothed_route = splev(u, tck)
    
    #Compute radius of curvature
    derivative1 = splev(u, tck, der=1)
    derivative2 = splev(u, tck, der=2) 
    xprim = derivative1[0]
    xprimprim = derivative2[0]
    yprim = derivative1[1]
    yprimprim = derivative2[1]
    curvature = 1.0 / np.abs((xprim*yprimprim - yprim* xprimprim) / np.power(xprim*xprim + yprim*yprim, 3 / 2))
    minCurvature = np.min(curvature)
    print("s is %g => Minimum curvature radius is %g"%(s,np.min(curvature)))

    #Perform bisection
    if minCurvature > minCurvatureConstraint:
        s_1 = s
        minCurvature_1 = minCurvature
    else:
        s_0 = s
        minCurvature_0 = minCurvature

Это может потребовать некоторых доработок, таких как итерации, чтобы найти подходящий s_1, но работает.