c++ - pierre - Comprobando si dos curvas de Bézier cúbicas se intersecan
graficar curvas de bezier online (8)
Para un proyecto personal, necesito averiguar si se intersecan dos curvas cúbicas de Bézier. No necesito saber dónde: solo necesito saber si lo hacen. Sin embargo, tendría que hacerlo rápido.
He estado limpiando el lugar y encontré varios recursos. Sobre todo, hay una pregunta aquí que tuvo una respuesta prometedora.
Entonces, después de averiguar qué es una matriz de Sylvester , qué es un determinant , qué es un resultant y por qué es útil , pensé que pensaba cómo funciona la solución. Sin embargo, la realidad difiere, y no funciona tan bien.
Jugando alrededor
He usado mi calculadora gráfica para dibujar dos splines Bézier (que llamaremos B 0 y B 1 ) que se intersecan. Sus coordenadas son las siguientes (P 0 , P 1 , P 2 , P 3 ):
(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)
El resultado es el siguiente, siendo B 0 la curva "horizontal" y B 1 la otra:
Siguiendo las instrucciones de la respuesta más votada de la pregunta mencionada, he restado B 0 a B 1 . Me dejó con dos ecuaciones (los ejes X e Y) que, según mi calculadora, son:
x = 9t^3 - 9t^2 - 3t + 2
y = 9t^3 - 9t^2 - 6t + 4
La matriz de Sylvester
Y a partir de eso he construido la siguiente matriz de Sylvester:
9 -9 -3 2 0 0
0 9 -9 -3 2 0
0 0 9 -9 -3 2
9 -9 -6 4 0 0
0 9 -9 -6 4 0
0 0 9 -9 -6 4
Después de eso, hice una función de C ++ para calcular los determinantes de matrices usando la expansión de Laplace :
template<int size>
float determinant(float* matrix)
{
float total = 0;
float sign = 1;
float temporaryMatrix[(size - 1) * (size - 1)];
for (int i = 0; i < size; i++)
{
if (matrix[i] != 0)
{
for (int j = 1; j < size; j++)
{
float* targetOffset = temporaryMatrix + (j - 1) * (size - 1);
float* sourceOffset = matrix + j * size;
int firstCopySize = i * sizeof *matrix;
int secondCopySize = (size - i - 1) * sizeof *matrix;
memcpy(targetOffset, sourceOffset, firstCopySize);
memcpy(targetOffset + i, sourceOffset + i + 1, secondCopySize);
}
float subdeterminant = determinant<size - 1>(temporaryMatrix);
total += matrix[i] * subdeterminant * sign;
}
sign *= -1;
}
return total;
}
template<>
float determinant<1>(float* matrix)
{
return matrix[0];
}
Parece funcionar bastante bien en matrices relativamente pequeñas (2x2, 3x3 y 4x4), así que espero que funcione también en matrices de 6x6. Sin embargo, no realicé pruebas exhaustivas y existe la posibilidad de que se rompa.
El problema
Si entendí correctamente la respuesta de la otra pregunta, el determinante debería ser 0 ya que las curvas se intersecan. Sin embargo, alimentando a mi programa la matriz de Sylvester que hice anteriormente, es -2916.
¿Es un error en mi final o en su final? ¿Cuál es la forma correcta de descubrir si se intersecan dos curvas cúbicas de Bézier?
¿Es un error en mi final o en su final?
¿Está basando su interpretación del determinante en el cuarto comentario adjunto a esta respuesta ? Si es así, creo que ahí radica el error. Reproduciendo el comentario aquí:
Si el determinante es cero, hay una raíz en X e Y en * exactamente el mismo valor de t, por lo que hay una intersección de las dos curvas. (aunque la t puede no estar en el intervalo 0..1).
No veo ningún problema con esta parte, pero el autor continúa diciendo:
Si el determinante es <> cero, puede estar seguro de que las curvas no se intersecan en ninguna parte.
No creo que eso sea correcto. Es perfectamente posible que las dos curvas se crucen en una ubicación donde los valores de t difieran, y en ese caso, habrá una intersección aunque la matriz tenga un determinante distinto de cero. Creo que esto es lo que está pasando en tu caso.
De ninguna manera soy un experto en este tipo de cosas, pero sigo un buen blog que habla mucho sobre curvas. Tiene un enlace a dos buenos artículos que hablan sobre su problema (el segundo enlace tiene una demostración interactiva y algún código fuente). Otras personas pueden tener una mejor comprensión del problema, pero espero que esto ayude.
Es este un problema difícil. Dividiría cada una de las 2 curvas de Bezier en, por ejemplo, 5-10 segmentos de línea discretos, y solo haría intersecciones de línea de línea.
foreach SampledLineSegment line1 in Bezier1
foreach SampledLineSegment line2 in Bezier2
if( line1 intersects line2 )
then Bezier1 intersects Bezier2
La intersección de las curvas de Bézier se realiza mediante el lenguaje de gráficos vectoriales de Asymptote (muy interesante): busque intersect() here .
Aunque no explican el algoritmo que realmente usan allí, excepto para decir que es de p. 137 de "The Metafont Book", parece que la clave de esto son dos propiedades importantes de las curvas de Bezier (que se explican en otras partes de ese sitio, aunque no puedo encontrar la página en este momento):
- Una curva Bezier siempre está contenida dentro del cuadro delimitador definido por sus 4 puntos de control
- Una curva Bezier siempre se puede subdividir en un valor t arbitrario en 2 curvas sub-Bezier
Con estas dos propiedades y un algoritmo para la intersección de polígonos, puede realizar una precisión arbitraria:
bezInt (B 1 , B 2 ):
- ¿Bbox (B 1 ) se cruza con Bbox (B 2 )?
- No: devuelve falso.
- Sí: continuar.
- ¿Es area (bbox (B 1 )) + area (bbox (B 2 )) <umbral?
- Sí: Devuelve verdadero.
- No: Continuar.
- Dividir B 1 en B 1a y B 1b en t = 0.5
- Dividir B 2 en B 2a y B 2b en t = 0.5
- Regresar a bezInt (B 1a , B 2a ) || bezInt (B 1a , B 2b ) || bezInt (B 1b , B 2a ) || bezInt (B 1b , B 2b ).
Esto será rápido si las curvas no se intersecan. ¿Es ese el caso habitual?
[EDITAR] Parece que el algoritmo para dividir una curva Bezier en dos se llama algoritmo de Casteljau .
No sé qué tan rápido será, pero si tiene dos curvas C1 (t) y C2 (k) se intersecan si C1 (t) == C2 (k). Entonces tienes dos ecuaciones (por x y por y) para dos variables (t, k). Puede resolver el sistema utilizando métodos numéricos con suficiente precisión. Cuando haya encontrado t, k parámetros, debe comprobar si hay un parámetro en [0, 1]. Si es que se intersecan en [0, 1].
Sí, lo sé, este hilo es aceptado y cerrado durante mucho tiempo, pero ...
Primero, me gustaría agradecerte, zneak , por una inspiración. Tu esfuerzo permite encontrar el camino correcto.
En segundo lugar, no estaba del todo contento con la solución aceptada. Este tipo se usa en mi freeware favorito IPE, y su bugzilla está lleno de quejas por una baja precisión y confiabilidad sobre un problema de intersección, entre ellos.
El truco faltante que permite resolver el problema de una manera propuesta por zneak : es suficiente para acortar una de las curvas en un factor k > 0, entonces el determinante de la matriz de Sylvester será igual a cero. Es obvio que si una curva acortada se intersecará, entonces la original también lo hará. Ahora, la tarea se cambia a la búsqueda de un valor adecuado para el factor k . Esto ha llevado al problema de resolver un polinomio univariado de 9 grados. Las raíces reales y positivas de este polinomio son responsables de los posibles puntos de intersección. (Esto no debería ser una sorpresa, dos curvas Bezier cúbicas pueden intersecar hasta 9 veces). La selección final se realiza para encontrar solo esos factores k , que proporcionan 0 <= t <= 1 para ambas curvas.
Ahora el código de Maxima, donde se establece el punto de inicio de 8 puntos proporcionado por zneak :
p0x:1; p0y:1;
p1x:2; p1y:4;
p2x:3; p2y:4;
p3x:4; p3y:3;
q0x:3; q0y:5;
q1x:3; q1y:6;
q2x:0; q2y:1;
q3x:3; q3y:1;
c0x:p0x;
c0y:p0y;
c1x:3*(p1x-p0x);
c1y:3*(p1y-p0y);
c2x:3*(p2x+p0x)-6*p1x;
c2y:3*(p2y+p0y)-6*p1y;
c3x:3*(p1x-p2x)+p3x-p0x;
c3y:3*(p1y-p2y)+p3y-p0y;
d0x:q0x;
d0y:q0y;
d1x:3*(q1x-q0x);
d1y:3*(q1y-q0y);
d2x:3*(q2x+q0x)-6*q1x;
d2y:3*(q2y+q0y)-6*q1y;
d3x:3*(q1x-q2x)+q3x-q0x;
d3y:3*(q1y-q2y)+q3y-q0y;
x:c0x-d0x + (c1x-d1x*k)*t+ (c2x-d2x*k^2)*t^2+ (c3x-d3x*k^3)*t^3;
y:c0y-d0y + (c1y-d1y*k)*t+ (c2y-d2y*k^2)*t^2+ (c3y-d3y*k^3)*t^3;
z:resultant(x,y,t);
En este punto, Maxima respondió:
(%o35)−2*(1004*k^9−5049*k^8+5940*k^7−1689*k^6+10584*k^5−8235*k^4−2307*k^3+1026*k^2+108*k+76)
Deja que Maxima resuelva esta ecuación:
rr: float( realroots(z,1e-20))
La respuesta es:
(%o36) [k=−0.40256438624399,k=0.43261490325108,k=0.84718739982868,k=2.643321910825066,k=2.71772491293651]
Ahora el código para seleccionar un valor correcto de k :
for item in rr do (
evk:ev(k,item),
if evk>0 then (
/*print("k=",evk),*/
xx:ev(x,item), rx:float( realroots(xx,1e-20)),/*print("x(t)=",xx," roots: ",rx),*/
yy:ev(y,item), ry:float( realroots(yy,1e-20)),/*print("y(t)=",yy," roots: ",ry),*/
for it1 in rx do ( t1:ev(t,it1),
for it2 in ry do ( t2:ev(t,it2),
dt:abs(t1-t2),
if dt<1e-10 then (
/*print("Common root=",t1," delta t=",dt),*/
if (t1>0) and (t1<=1) then ( t2:t1*evk,
if (t2>0) and (t2<=1) then (
x1:c0x + c1x*t1+ c2x*t1^2+ c3x*t1^3,
y1:c0y + c1y*t1+ c2y*t1^2+ c3y*t1^3,
print("Intersection point: x=",x1, " y=",y1)
)))))/*,disp ("-----")*/
));
La respuesta de Maxima:
"Intersection point: x="1.693201254437358" y="2.62375005067273
(%o37) done
Aunque no solo hay una miel. El método presentado es inutilizable si:
- P0 = Q0, o más generalmente, si P0 se encuentra en la segunda curva (o su extensión). Uno puede intentar cambiar las curvas.
- casos extremadamente raros, cuando ambas curvas pertenecen a una familia K (por ejemplo, sus extensiones infinitas son las mismas).
- mantener un ojo en (sqr (c3x) + sqr (c3y)) = 0 o (sqr (d3x) + sqr (3y)) = 0 casos, aquí una cuadrática pretende ser una curva Bezier cúbica.
Uno puede preguntar, ¿por qué un acortamiento se realiza una sola vez. Es suficiente debido a la ley inversa inversa, que se descubrió de forma pasiva , pero esta es otra historia.
Si está haciendo esto para el código de producción, sugeriría el algoritmo de recorte Bezier. Está bien explicado en la cagd.cs.byu.edu/~557/text/ch7.pdf (pdf), funciona para cualquier grado de curva de Bezier y es rápido y robusto.
Si bien el uso de buscadores de raíz o matrices estándar puede ser más directo desde un punto de vista matemático, el recorte de Bezier es relativamente fácil de implementar y depurar, y en realidad tendrá menos error de punto flotante. Esto se debe a que cada vez que se crean nuevos números, se realizan promedios ponderados (combinaciones convexas), por lo que no hay posibilidad de extrapolar según las entradas ruidosas.
Yo diría que la respuesta más fácil y rápida es subdividirlas en líneas muy pequeñas y encontrar los puntos donde se intersecan las curvas, si es que realmente lo hacen.
public static void towardsCubic(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double t) {
double x, y;
x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;
xy[0] = x;
xy[1] = y;
}
public static void towardsQuad(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double t) {
double x, y;
x = (1 - t) * (1 - t) * x0 + 2 * (1 - t) * t * x1 + t * t * x2;
y = (1 - t) * (1 - t) * y0 + 2 * (1 - t) * t * y1 + t * t * y2;
xy[0] = x;
xy[1] = y;
}
Obviamente, la respuesta de fuerza bruta es mala. Vea la respuesta de bo {4}, y hay mucha geometría lineal y detección de colisiones que realmente ayudarán bastante.
Elija el número de cortes que desee para las curvas. Algo como 100 debería ser genial.
Tomamos todos los segmentos y los ordenamos en función del mayor valor de y que tienen. Luego, agregamos una segunda bandera en la lista para el valor más pequeño de y para ese segmento.
Mantenemos una lista de bordes activos.
Recorremos la lista de segmentos ordenada por y, cuando encontramos un segmento inicial, lo agregamos a la lista activa. Cuando alcanzamos el valor marcado con una pequeña y, eliminamos ese segmento de la lista activa.
Luego, podemos simplemente iterar a través de todo el conjunto de segmentos con lo que equivale a una línea de escaneo, aumentando nuestra y monótonamente a medida que simplemente iteramos la lista. Recorremos los valores en nuestra lista ordenada, que normalmente eliminará un segmento y agregará un nuevo segmento (o para dividir y fusionar nodos, agregar dos segmentos o eliminar dos segmentos). Manteniendo una lista activa de segmentos relevantes.
Ejecutamos la comprobación de intersección de falla rápida a medida que agregamos un nuevo segmento activo a la lista de segmentos activos, solo contra ese segmento y los segmentos actualmente activos.
Por lo tanto, en todo momento sabemos exactamente qué segmentos de línea son relevantes, a medida que recorremos los segmentos muestreados de nuestras curvas. Sabemos que estos segmentos comparten la superposición en los y-coords. Podemos fallar rápidamente cualquier segmento nuevo que no se superponga con respecto a sus x-coords. En el raro caso de que se superpongan en las x-coords, verificamos si estos segmentos se intersecan.
Esto probablemente reducirá el número de verificaciones de intersecciones de línea básicamente a la cantidad de intersecciones.
foreach(segment in sortedSegmentList) {
if (segment.isLeading()) {
checkAgainstActives(segment);
actives.add(segment);
}
else actives.remove(segment)
}
checkAgainstActive () simplemente verificará si este segmento y cualquier segmento en la lista activa se superponen con sus x-coords, si lo hacen, luego ejecutan la línea en el cruce de intersección en ellos, y toman la acción apropiada.
También tenga en cuenta que esto funcionará para cualquier número de curvas, cualquier tipo de curvas, cualquier orden de curvas, en cualquier mezcla. Y si recorremos la lista completa de segmentos, encontrará todas las intersecciones. Encontrará todos los puntos en los que Bezier interseca un círculo o cada intersección que una docena de curvas Bezier cuadráticas tienen entre sí (o ellas mismas), y todas en la misma fracción de segundo.
Las notas del cagd.cs.byu.edu/~557/text/ch7.pdf menudo citadas, para el algoritmo de subdivisión:
"Una vez que un par de curvas se ha subdividido lo suficiente como para que cada una pueda aproximarse por un segmento de línea a una tolerancia"
Podemos, literalmente, simplemente saltarse al intermediario. Podemos hacer esto lo suficientemente rápido para simplemente comparar segmentos de línea con una cantidad tolerable de error de la curva. Al final, eso es lo que hace la respuesta típica.
En segundo lugar, tenga en cuenta que la mayor parte del aumento de velocidad para la detección de colisión aquí, es decir, la lista ordenada de segmentos ordenados en función de su y más alta para agregar a la lista activa, y la y más baja para eliminar de la lista activa, también se puede hacer Para los polígonos del casco de la curva Bezier directamente. Nuestro segmento de línea es simplemente un polígono de orden 2, pero podríamos hacer cualquier número de puntos de manera trivial, y marcar el cuadro delimitador de todos los puntos de control en el orden de curva que estemos tratando. Entonces, en lugar de una lista de segmentos activos, tendríamos una lista de puntos de polígonos de casco activos. Podríamos simplemente usar el algoritmo de De Casteljau para dividir las curvas, y así muestrearlas como curvas subdivididas en lugar de segmentos de línea. Entonces, en lugar de 2 puntos, tendríamos 4 o 7 o lo que sea, y ejecutaremos la misma rutina siendo bastante aptos para fallar rápidamente.
Agregando el grupo relevante de puntos al máximo y, eliminándolo al mínimo y comparando solo la lista activa. Por lo tanto, podemos implementar mejor y rápidamente el algoritmo de subdivisión Bezier. Simplemente encontrando la superposición del cuadro delimitador, y luego subdividiendo las curvas que se superponen, y eliminando las que no lo hicieron. Como, de nuevo, podemos hacer cualquier número de curvas, incluso las subdivididas de las curvas en la iteración anterior. Y al igual que nuestra aproximación de segmento, resuelva rápidamente para cada posición de intersección entre cientos de curvas diferentes (incluso de órdenes diferentes) muy rápidamente. Simplemente revise todas las curvas para ver si los cuadros delimitadores se superponen, si lo hacen, los subdividimos, hasta que nuestras curvas sean lo suficientemente pequeñas o salimos corriendo de ellas.