Evaluación de Bezier en O (n) -- java campo con performance campo con algorithm camp codereview Relacionados El problema

Bezier evaluation in O(n)


4
vote

problema

Español

La forma realmente estúpida de evaluar las curvas de Bezier es con la recursión, que es $ O (2 ^ N) $, que se puede bajar a $ O (n ^ 2) $ con memoización. El algoritmo de De Casteljau mejora esto, y sigue siendo todavía $ O (n ^ 2) $, pero más rápido.

Mal interpreté a De Casteljau y pensé que era la forma de Bernstein:

$$ sum_ {i = 0} ^ {n-1} { binom {n} {i} p_i (1-t) ^ {n-i-1} T ^ i} $$

Donde hay puntos de control $ N $, incluido el inicio y finalizar en nombre $ P_0, P_1, P_2, DOTS, P_ {N-1} $.

Cálculo de la combinación es normalmente $ O (n) $, que haría este algoritmo en su totalidad $ O (n ^ 2) $, ya que hay $ n $ términos a suma, sin embargo, hay un truco:

$$ binom {n} {k} = prod_ {i = 1} ^ {k} { frac {n + 1-i} {k}} $$

Podemos almacenar el resultado después de cada iteración, lo que resulte en toda la secuencia calculada en el tiempo $ O (N) $, lo que hace que el algoritmo sea un completo $ O (n) $.

$ 2 leq n leq 5 $ es el caso de uso general (2 y 4 probablemente en realidad), por lo que está codificado. Para $ N GEQ 6 $, que utiliza este algoritmo, he hecho un número ridículo de optimizaciones.

El incremento 99887766655544330 es $ 2.2E {-16} $, y en función de algunos puntos de referencia rápidos con datos aleatorios, nunca fue más inexacta que $ 1E {- 14} $, por lo que la inexactitud es aceptable.

Por muy grande $ n $ (cientos, 1200 fue donde lo vi por primera vez) la función de combinación excedió el valor máximo para double NaN < / Código>, mientras que de Casteljau, aunque mucho más lento, devuelve lo que era más o menos el valor correcto.

El código final:

  public static double bezier2(double a,double b,double t){     //Total of 3 floating point operations     //      1 2  3     return a+t*(b-a); }  public static double bezier(double t,double... ds){     return bezier(ds,t); }  public static double bezier(double[] ds,double t){     int count = ds.length;     switch(count){     case 0:throw new IllegalArgumentException("Must have at least two items to interpolate between");     case 1:return ds[0];     case 2:return bezier2(ds[0],ds[1],t);     case 3:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double t1 = 1d - t;         /*          * Hardcoded for n=3          * Total of 8+1=9 floating point operations          *       1  2  3 4  5  6 7 8          */         return (a*t1+2d*b*t)*t1+c*t*t;     }     case 4:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double d = ds[3];         double t1 = 1d - t;         /*          * Hardcoded for n=4          * Total of 13+1=9 floating point operations          *       1  2  3 4  5  6  7   8 9 10 11 12 13          */         return (a*t1+3d*b*t)*t1*t1+(3d*c*t1+d*t)*t*t;     }     case 5:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double d = ds[3];         double e = ds[4];         double t1 = 1d - t;         double i = t1*t1;         double j = t*t;         double k = t1*t;         double l = 4d*k;         /*          * Hardcoded for n=5          * Total of 12+5=17 floating point operations          *       1  2  3  4   5 6   7   8   9 10 11  12          */         return (a*i + b*l + 6d*c*j) * i + (d*l + e*j) * j;     }     case 6:     case 7:     case 8:     case 9:     case 10:     case 11:return decasteljauBezier(ds,t);     default:{         double t1 = 1d - t;         int n1 = count - 1;         int halfn = (n1>>1)+1;         int halfn1 = halfn+1;         double[] choose;         if(count<29){             choose = new double[halfn1];             int[] chooseInt = chooseIntRange(n1,halfn);             for(int i=0;i<halfn1;i++){                 choose[i]=chooseInt[i];             }         }else if(count<60){             choose = new double[halfn1];             long[] chooseLong = chooseLongRange(n1,halfn);             for(int i=0;i<halfn1;i++){                 choose[i]=chooseLong[i];             }         }else{             choose = chooseDoubleRange(n1,halfn);         }         double[] terms = new double[count];         double power = 1d;         for(int i=0;i<halfn;i++){             terms[i] = ds[i] * choose[i] * power;             power *= t;         }         for(int i=halfn;i<count;i++){             terms[i] = ds[i] * choose[n1-i] * power;             power *= t;         }         power = t1;         for(int i=1;i<count;i++){             terms[n1-i] *= power;             power *= t1;         }         double sum = 0d;         for(double v:terms){             sum += v;         }         return sum;     }     } }  public static double decasteljauBezier(double[] ds,double t){     int n = ds.length;     double[] result = Arrays.copyOf(ds, n);     for(int i=n-1;i>0;i--){         for(int j=0;j<i;j++){             result[j]+=(result[j+1]-result[j])*t;         }     }     return result[0]; }  public static int[] chooseIntRange(int n,int k){     int n1 = n+1;     int product = 1;     int[] result = new int[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = product*(n1-i)/i;         result[i] = product;     }     return result; }  public static long[] chooseLongRange(int n,int k){     int n1 = n+1;     int product = 1;     long[] result = new long[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = product*(n1-i)/i;         result[i] = product;     }     return result; }  public static double[] chooseDoubleRange(int n,int k){     int n1 = n+1;     double product = 1d;     double[] result = new double[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = (product * (n1-i)) / i;         result[i] = product;     }     return result; }   

Notas:

  1. Esto fue optimizado a mano por mí, y probablemente nuevamente por el compilador. Explicaría cómo $ 2 leq n leq 5 $ fue de aproximadamente la misma velocidad, pero $ n = 6 $, de repente, fue $ frac {1} {3} $ de la velocidad.

  2. Estos números no se seleccionan aleatoriamente. Se basan en los valores máximos para int S y 9988777665544335 s. Puede ser posible aumentar estos un bit (y, por lo tanto, mejorar el algoritmo que sea tan ligeramente), pero me gusta la seguridad.

  3. Semi-in-Place de Casteljau's con todas las optimizaciones que podría pensar. Es bastante rápido pero sigue siendo $ O (n ^ 2) $.

de Casteljau fue más rápido por $ 6 leq n leq 11 $, así que pensé en redirigir. Pero como resulta, el costo adicional de esa declaración única 9988776655544336 y una llamada de redireccionamiento significaba que solo valía la pena redirigir por $ 6 leq n leq 9 $. Incluso esa redirección podría desaparecer en el futuro con más optimización. No es como si el tiempo no fuera comparable por un pequeño "$ n $.

Beating $ O (n) $ es definitivamente imposible, pero estoy seguro de que hay un mejor algoritmo, o si no, una mejor manera de implementar este. Si no, ni siquiera eso, seguramente hay optimizaciones que he perdido.

Original en ingles

The really stupid way to evaluate bezier curves is with recursion, which is \$O(2^n)\$, which can be lowered to \$O(n^2)\$ with memoization. De Casteljau's algorithm improves on this, and is still \$O(n^2)\$ but faster.

I misinterpreted De Casteljau's and thought it was the Bernstein form:

$$\sum_{i=0}^{n-1}{\binom{n}{i}p_i (1-t)^{n-i-1} t^i}$$

Where there are \$n\$ control points including the start and end named \$p_0, p_1,p_2, \dots, p_{n-1}\$.

Calculating the combination is normally \$O(n)\$ which would make this algorithm as a whole \$O(n^2)\$ since there are \$n\$ terms to sum, however, there is a cheat:

$$\binom{n}{k}=\prod_{i=1}^{k}{\frac{n+1-i}{k}}$$

We can store the result after each iteration, resulting in the entire sequence calculated in \$O(n)\$ time, making the algorithm as a whole \$O(n)\$.

\$2\leq n\leq 5\$ is the general use case (2 and 4 probably actually), so it's hardcoded. For \$n\geq 6\$, which uses this algorithm, I've done a ridiculous number of optimizations.

The smallest double increment (one mantissa bit) is \$2.2E{-16}\$, and based on some quick benchmarks with random data it never was more inaccurate than \$1E{-14}\$, so the inaccuracy is acceptable.

For very large \$n\$ (hundreds, 1200 was where I first saw it) the combination function exceeded the max value for doubles, causing the algorithm to return NaN, whereas De Casteljau's, though much slower, returned what was more or less the correct value.

The final code:

public static double bezier2(double a,double b,double t){     //Total of 3 floating point operations     //      1 2  3     return a+t*(b-a); }  public static double bezier(double t,double... ds){     return bezier(ds,t); }  public static double bezier(double[] ds,double t){     int count = ds.length;     switch(count){     case 0:throw new IllegalArgumentException("Must have at least two items to interpolate between");     case 1:return ds[0];     case 2:return bezier2(ds[0],ds[1],t);     case 3:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double t1 = 1d - t;         /*          * Hardcoded for n=3          * Total of 8+1=9 floating point operations          *       1  2  3 4  5  6 7 8          */         return (a*t1+2d*b*t)*t1+c*t*t;     }     case 4:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double d = ds[3];         double t1 = 1d - t;         /*          * Hardcoded for n=4          * Total of 13+1=9 floating point operations          *       1  2  3 4  5  6  7   8 9 10 11 12 13          */         return (a*t1+3d*b*t)*t1*t1+(3d*c*t1+d*t)*t*t;     }     case 5:{         double a = ds[0];         double b = ds[1];         double c = ds[2];         double d = ds[3];         double e = ds[4];         double t1 = 1d - t;         double i = t1*t1;         double j = t*t;         double k = t1*t;         double l = 4d*k;         /*          * Hardcoded for n=5          * Total of 12+5=17 floating point operations          *       1  2  3  4   5 6   7   8   9 10 11  12          */         return (a*i + b*l + 6d*c*j) * i + (d*l + e*j) * j;     }     case 6:     case 7:     case 8:     case 9:     case 10:     case 11:return decasteljauBezier(ds,t);     default:{         double t1 = 1d - t;         int n1 = count - 1;         int halfn = (n1>>1)+1;         int halfn1 = halfn+1;         double[] choose;         if(count<29){             choose = new double[halfn1];             int[] chooseInt = chooseIntRange(n1,halfn);             for(int i=0;i<halfn1;i++){                 choose[i]=chooseInt[i];             }         }else if(count<60){             choose = new double[halfn1];             long[] chooseLong = chooseLongRange(n1,halfn);             for(int i=0;i<halfn1;i++){                 choose[i]=chooseLong[i];             }         }else{             choose = chooseDoubleRange(n1,halfn);         }         double[] terms = new double[count];         double power = 1d;         for(int i=0;i<halfn;i++){             terms[i] = ds[i] * choose[i] * power;             power *= t;         }         for(int i=halfn;i<count;i++){             terms[i] = ds[i] * choose[n1-i] * power;             power *= t;         }         power = t1;         for(int i=1;i<count;i++){             terms[n1-i] *= power;             power *= t1;         }         double sum = 0d;         for(double v:terms){             sum += v;         }         return sum;     }     } }  public static double decasteljauBezier(double[] ds,double t){     int n = ds.length;     double[] result = Arrays.copyOf(ds, n);     for(int i=n-1;i>0;i--){         for(int j=0;j<i;j++){             result[j]+=(result[j+1]-result[j])*t;         }     }     return result[0]; }  public static int[] chooseIntRange(int n,int k){     int n1 = n+1;     int product = 1;     int[] result = new int[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = product*(n1-i)/i;         result[i] = product;     }     return result; }  public static long[] chooseLongRange(int n,int k){     int n1 = n+1;     int product = 1;     long[] result = new long[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = product*(n1-i)/i;         result[i] = product;     }     return result; }  public static double[] chooseDoubleRange(int n,int k){     int n1 = n+1;     double product = 1d;     double[] result = new double[k+1];     result[0] = 1;     for(int i=1;i<=k;i++){         product = (product * (n1-i)) / i;         result[i] = product;     }     return result; } 

Notes:

  1. This was optimized by hand by me, and likely again by the compiler. It would explain how \$2\leq n\leq 5\$ was all roughly the same speed but \$n=6\$ then suddenly was \$\frac{1}{3}\$ of the speed.

  2. These numbers aren't randomly chosen. They're based on the max values for ints and longs. It might be possible to increase these a bit (and therefore improve the algorithm ever so slightly) but I like safety.

  3. Semi-in-place De Casteljau's with all the optimizations I could think of. It's pretty fast but still \$O(n^2)\$.

De Casteljau's was faster for \$6\leq n\leq 11\$, so I thought of redirecting. But as it turns out, the extra cost of that single if statement and a redirect call meant that it was only worth redirecting for \$6\leq n\leq 9\$. Even that redirect might be gone in the future with more optimizing. It's not like the time isn't comparable for such a small \$n\$.

Beating \$O(n)\$ is definitely impossible, but I'm certain there is a better algorithm out there, or if not, a better way to implement this one. If not even that, surely there are optimizations I've missed.

        

Lista de respuestas

2
 
vote
vote
La mejor respuesta
 

Debe mirar sus pruebas de unidad, porque hay un error malo en chooseLongRange . int product = 1 debe ser long product = 1 .


      if(count<29){//See note 1         choose = new double[halfn1];         int[] chooseInt = chooseIntRange(n1,halfn);         for(int i=0;i<halfn1;i++){             choose[i]=chooseInt[i];         }     }else if(count<60){//See note 2         choose = new double[halfn1];         long[] chooseLong = chooseLongRange(n1,halfn);         for(int i=0;i<halfn1;i++){             choose[i]=chooseLong[i];         }     }else{         choose = chooseDoubleRange(n1,halfn);     }   

imo Esto es excesivamente complicado. A menos que esté reutilizando chooseIntRange chooseLongRange5 en contextos que requieren int[] o long[] , ¿por qué no hacer que vuelvan < Código> double[] ? Y luego puede presionar la caja dividida para que se use el tipo, exponiendo un solo método y dejando la responsabilidad en el método en el que realmente pertenece.


Estoy ligeramente desconcertado por ", pero según lo resulta, el costo adicional de esa sola declaración de SI SIGNO y una llamada de redireccionamiento significaba que solo valía la pena redirigir por 6≤n≤9 ". ¿Qué "Single IF Declaración"? ¿Te refieres a la sobrecarga del switch Búsqueda de la tabla?


Tengo la impresión de que le preocupa más sobre el rendimiento que el análisis numérico, pero creo que vale la pena hacer este punto de todos modos.

Hay una diferencia importante entre

  int product = 10  

y

  int product = 11  

Si sustituye $ (a, b, t) = (b, a, 1-t) $, entonces el segundo le dará el mismo resultado, pero la primera no necesariamente. En particular, cuando $ T aprox 1 $ y $ b ll debería estar muy cerca de $ b $, pero en la primera expresión la subexpression $ b - a $ perderá La precisión de $ B $ y no se puede recuperar en el resultado final.

 

You need to look at your unit tests, because there's a bad bug in chooseLongRange. int product = 1 should be long product = 1.


    if(count<29){//See note 1         choose = new double[halfn1];         int[] chooseInt = chooseIntRange(n1,halfn);         for(int i=0;i<halfn1;i++){             choose[i]=chooseInt[i];         }     }else if(count<60){//See note 2         choose = new double[halfn1];         long[] chooseLong = chooseLongRange(n1,halfn);         for(int i=0;i<halfn1;i++){             choose[i]=chooseLong[i];         }     }else{         choose = chooseDoubleRange(n1,halfn);     } 

IMO this is excessively complicated. Unless you're reusing chooseIntRange and chooseLongRange in contexts which require int[] or long[], why not make them return double[]? And then you can push in the case split for the type to use, exposing a single method and leaving the responsibility in the method where it really belongs.


I'm slightly puzzled by "But as it turns out, the extra cost of that single if statement and a redirect call meant that it was only worth redirecting for 6xe2x89xa4nxe2x89xa49". What "single if statement"? Do you mean the overhead of the switch table lookup?


I get the impression that you care more about performance than numerical analysis, but I think it's worth making this point anyway.

There's an important difference between

a + (b - a) * t 

and

a * (1 - t) + b * t 

If you substitute \$(a, b, t) = (b, a, 1-t)\$ then the second will give you the same result, but the first won't necessarily. In particular, when \$t \approx 1\$ and \$b \ll a\$ the result should be very close to \$b\$, but in the first expression the subexpression \$b - a\$ will lose the precision of \$b\$ and it can't be recovered in the final result.

 
 

Relacionados problema

35  Demasiados bucles en la aplicación de dibujo  ( Too many loops in drawing app ) 
Tengo un método que tiene muchos bucles: #ifndef __RUNES_STRUCTURES_H #define __RUNES_STRUCTURES_H /* Runes structures. */ struct Game { char board[2...

25  Algoritmo para transformar una palabra a otra a través de palabras válidas  ( Algorithm to transform one word to another through valid words ) 
He estado practicando retroceso y quería saber cómo puedo mejorar mi código. Por ejemplo, no quiero usarlo global. Además, no estoy seguro de si mi código fun...

5  Orden de número más grande en cadena  ( Largest number order in string ) 
Dada una cadena, suponiendo que la cadena sea solo números, reorganice la cadena a la que sea el mayor número posible. a continuación es mi solución al pr...

56  Proyecto Euler Problema 1 en Python - Múltiples de 3 y 5  ( Project euler problem 1 in python multiples of 3 and 5 ) 
Me gustaría sugerencias para optimizar esta solución de fuerza bruta a problema 1 . El algoritmo actualmente comprueba cada entero entre 3 y 1000. Me gustarí...

2  Dos formas de aleatorias aleatoriamente las tarjetas  ( Two ways to randomly shuffle cards ) 
Aquí hay dos implementaciones que escribí para aleatorizar las tarjetas. El primer método ( dt5 ) Selecciona una tarjeta aleatoria, luego lo quita al frent...

1  Retire todos los nodos que no se encuentren en ningún camino con suma> = k  ( Remove all nodes which dont lie in any path with sum k ) 
Dado un árbol binario, una ruta completa se define como un camino desde la raíz a una hoja. La suma de todos los nodos en ese camino se define como la suma d...

6  Encontrar el siguiente palíndromo de una cadena de números  ( Finding the next palindrome of a number string ) 
Aquí está el problema: Un entero positivo se llama palíndromo si su representación en el El sistema decimal es el mismo cuando se lee de izquierda a dere...

8  Simple GCD Utility en Java  ( Simple gcd utility in java ) 
i anteriormente discutido El rendimiento se refiere a diferentes algoritmos GCD. Escribí una simple clase de Java que implementa el algoritmo binario GCD. E...

1  Compruebe si dos cadenas son permutación entre sí  ( Check if two strings are permutation of each other ) 
private String sort(String word) { char[] content = word.toCharArray(); Arrays.sort(content); return new String(content); } private boolea...

5  Proyecto EULER NO. 17: contando letras para escribir los números de 1 a 1000  ( Project euler no 17 counting letters to write the numbers from 1 to 1000 ) 
Soy muy nuevo en la programación y, cierto, estoy avergonzado de compartir mi código para la crítica. Este código funciona y produce la respuesta correcta a l...




© 2022 respuesta.top Reservados todos los derechos. Centro de preguntas y respuestas reservados todos los derechos