Datos aproximados con una curva bezier cúbica multisegmento y una distancia, así como una contraposición de curvatura

Tengo algunos datos geográficos (la imagen de abajo muestra la trayectoria de un río como puntos rojos), que quiero aproximar utilizando una curva bezier cúbica de múltiples segmentos. A través de otras preguntas sobre stackoverflow aquí y aquí , encontré el algoritmo de Philip J. Schneider de “Graphics Gems”. Lo implementé con éxito y puedo informar que incluso con miles de puntos es muy rápido. Desafortunadamente, la velocidad tiene algunas desventajas, a saber, que la adaptación se hace bastante descuidadamente. Considere el siguiente gráfico:

curva de bezier multi segmento

Los puntos rojos son mis datos originales y la línea azul es el segmento multiseguro más pequeño creado por el algoritmo de Schneider. Como puede ver, la entrada al algoritmo fue una tolerancia que es al menos tan alta como indica la línea verde. Sin embargo, el algoritmo crea una curva de bezier que tiene demasiados giros bruscos. También ves estos innecesarios giros bruscos en la imagen. Es fácil imaginar una curva Bezier con giros menos nítidos para los datos mostrados mientras se mantiene la condición de tolerancia máxima (solo presione la curva Bezier un poco en la dirección de las flechas magenta). El problema parece ser que el algoritmo selecciona puntos de datos de mis datos originales como puntos finales de las curvas de bezier individuales (los puntos de flechas magenta indican algunos sospechosos). Con los puntos finales de las curvas de Bézier restringidos de esa manera, está claro que el algoritmo a veces producirá curvaturas bastante pronunciadas.

Lo que estoy buscando es un algoritmo que se aproxime a mis datos con una curva más amplia de múltiples segmentos con dos restricciones:

  • la curva bezier de múltiples segmentos nunca debe estar a más de una cierta distancia de los puntos de datos (esto lo proporciona el algoritmo de Schneider)
  • la curva bezier multisegmento nunca debe crear curvaturas demasiado definidas. Una forma de verificar este criterio sería hacer un círculo con el radio de curvatura mínimo a lo largo de la curva de Bizier multisegmento y verificar si toca todas las partes de la curva a lo largo de su trayectoria. Aunque parece que hay un método mejor que involucra el producto cruzado de la primera y la segunda derivada

Las soluciones que encontré que crean mejores ajustes por desgracia funcionan solo para curvas bezier individuales (y omiten la pregunta de cómo encontrar buenos puntos de inicio y final para cada curva bezier en la curva bezier de segmentos múltiples) o no permiten un contraint de curvatura mínima. Siento que el contraint mínimo de curvatura es la condición difícil aquí.

Aquí otro ejemplo (esto es dibujado a mano y no es 100% preciso):

algunos ejemplos

Supongamos que la figura uno muestra tanto la restricción de curvatura (el círculo debe ajustarse a lo largo de toda la curva) como la distancia máxima de cualquier punto de datos desde la curva (que es el radio del círculo en verde). Una aproximación exitosa del camino rojo en la figura dos se muestra en azul. Esa aproximación respeta la condición de curvatura (el círculo puede rodar dentro de toda la curva y la toca en todas partes), así como la condición de distancia (que se muestra en verde). La figura tres muestra una aproximación diferente al camino. Si bien respeta la condición de distancia, está claro que el círculo ya no encaja en la curvatura. La figura cuatro muestra un camino que es imposible de aproximar con las restricciones dadas porque es demasiado puntiagudo. Este ejemplo debe ilustrar que para aproximar correctamente algunos giros puntiagudos en la ruta, es necesario que el algoritmo elija puntos de control que no forman parte de la ruta. La figura tres muestra que si se eligieron puntos de control a lo largo del camino, la restricción de curvatura ya no se puede cumplir. Este ejemplo también muestra que el algoritmo debe cerrarse en algunas entradas ya que no es posible aproximarlo con las restricciones dadas.

¿Existe una solución a este problema? La solución no tiene que ser rápida. Si se tarda un día en procesar 1000 puntos, está bien. La solución tampoco tiene que ser óptima en el sentido de que debe dar como resultado un ajuste de mínimos cuadrados.

Al final, implementaré esto en C y Python, pero también puedo leer la mayoría de los otros idiomas.

Encontré la solución que cumple con mi criterio. La solución es primero encontrar una B-Spline que se aproxime a los puntos en el sentido de mínimos cuadrados y luego convertir esa spline en una curva Bezier de múltiples segmentos. Las B-Splines tienen la ventaja de que, en contraste con las curvas de Bézier, no pasan a través de los puntos de control y proporcionan 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 python. Supongamos que leo mis datos en las listas 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 se ve así:

introduzca la descripción de la imagen aquí

Si quiero que la curva sea más suave, entonces puedo boost el parámetro s a splprep . Si quiero una aproximación más cercana a los datos, puedo disminuir el parámetro s para una menor suavidad. Al pasar por varios parámetros s 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 más pequeña. La respuesta en este correo electrónico por Zachary Pincus. Replicaré 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(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) 

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

  1. datos poligonales

    encuentre el orden de los puntos para que encuentre los puntos más cercanos entre sí e intente conectarlos ‘por líneas’. Evitar 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 se alcanza el mínimo o el 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 usa estos puntos como puntos de control. Recomiendo encarecidamente el polinomio de interpolación tanto para x como para 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 

    donde 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 usa 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 la curva de p1 a p2

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

[edit1] los bordes demasiado afilados son 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 … elegir qué es mejor para ti

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

    dx/dl o dy/dl donde 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 conduce a un radio de curvatura demasiado pequeño

    calcular la intersección de los segmentos de línea adyacentes (líneas negras) punto medio. Los ejes perpendiculares, como en la imagen (líneas rojas), la distancia entre ellos y el punto de unión (línea azul) es su radio de curvatura. Cuando el radio de curvatura es menor, entonces su límite elimina ese punto …

    radio de curvatura

    Ahora, si realmente solo necesitas los cúbicos BEZIER, puedes convertir mi interpolación cúbica a cúbica BEZIER así:

 // --------------------------------------------------------------------------- // 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; 

En caso de que necesite la conversión inversa ver:

  • Curva de Bezier con puntos de control dentro de la curva.