Aproximación de datos con una curva de Bézier cúbica de varios segmentos y una distancia, así como una restricción de curvatura

Solución:

Encontré la solución que cumple con mi criterea. La solución es encontrar primero un B-Spline que se aproxime a los puntos en el sentido de mínimos cuadrados y luego convertir ese spline en una curva Bézier de múltiples segmentos. Las B-Splines tienen la ventaja de que, a diferencia de las curvas Bézier, no pasarán por los puntos de control, además de proporcionar una forma de especificar la “suavidad” deseada de la curva de aproximación. La funcionalidad necesaria para generar tal spline se implementa en la biblioteca FITPACK a la que scipy ofrece un enlace de Python. Supongamos que leo mis datos en las listas x y y, entonces puedo hacer:

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()

El resultado entonces se ve así:

Si quiero que la curva sea más suave, puedo aumentar la s parámetro a splprep. Si quiero una aproximación más cercana a los datos, puedo disminuir la s parámetro para menos suavidad. Pasando por múltiples s parámetros programáticamente puedo encontrar un buen parámetro que se ajuste a los requisitos dados.

Sin embargo, la pregunta es cómo convertir ese resultado en una curva Bézier. La respuesta en este correo electrónico de Zachary Pincus. Reproduciré su solución aquí para dar una respuesta completa a mi pregunta:

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
  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)

Así que B-Splines, FITPACK, numpy y scipy me salvaron el día 🙂

1. poligonizar datos

   encuentra el orden de los puntos para que encuentres los puntos más cercanos entre sí y pruébalos para conectarlos ‘por líneas’. Evite volver al punto de origen

2. calcular la derivación a lo largo de la ruta

   es el cambio de dirección de las ‘líneas’ donde golpea el mínimo o máximo local allí está su punto de control … Haga esto para reducir sus datos de entrada (deje solo los puntos de control).

3. curva

   ahora use estos puntos como puntos de control. Recomiendo encarecidamente el polinomio de interpolación para ambos x y y por separado, por ejemplo, algo como esto:

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

   dónde a0..a3 se calculan así:

   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 se calculan de la misma manera pero usan coordenadas y, por supuesto
   • p0..p3 son puntos de control para la curva de interpolación cúbica
   • t =<0.0,1.0> es el parámetro de curva de p1 para p2

   esto asegura que la posición y la primera derivación sean continuas (c1) y también puede usar BEZIER, pero no será una coincidencia tan buena como esta.

[edit1] bordes demasiado afilados es un GRAN problema

Para resolverlo, puede eliminar puntos de su conjunto de datos antes de obtener los puntos de control. Puedo pensar en dos formas de hacerlo ahora mismo … elige lo que sea mejor para ti

1. eliminar puntos del conjunto de datos con una primera derivación demasiado alta

   dx/dl o dy/dl dónde x,y son coordenadas y l es la longitud de la curva (a lo largo de su trayectoria). El cálculo exacto del radio de curvatura a partir de la derivación de la curva es complicado

2. eliminar puntos del conjunto de datos que conducen a un radio de curvatura demasiado pequeño

   calcular la intersección de los segmentos de línea vecinos (líneas negras) en el punto medio. Ejes perpendiculares como en la imagen (líneas rojas) la distancia de la misma y el punto de unión (línea azul) es su radio de curvatura. Cuando el radio de curvatura es más pequeño, su límite elimina ese punto …

   

   ahora, si realmente necesita solo BEZIER cúbicos, puede convertir mi interpolación cúbica a BEZIER cúbica de esta manera:

//  ---------------------------------------------------------------------------
//  x=cx[0]+(t*cx[1])+(tt*cx[2])+(ttt*cx[3]); // cubic x=f
//  ---------------------------------------------------------------------------
//  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;

En caso de que necesite la conversión inversa, consulte:

• Curva de Bézier con puntos de control dentro de la curva.

---

---

Utiliza Nuestro Buscador