performance sorting math geometry angle

performance - La forma más rápida de ordenar los vectores por ángulo sin realmente calcular ese ángulo



sorting math (7)

Comencé a jugar un poco con esto y me di cuenta de que la especificación estaba incompleta. atan2 tiene una discontinuidad, porque como dx y dy varían, hay un punto en el que atan2 saltará entre -pi y + pi. La siguiente gráfica muestra las dos fórmulas sugeridas por @MvG, y de hecho ambas tienen la discontinuidad en un lugar diferente en comparación con atan2 . (NB: agregué 3 a la primera fórmula y 4 a la alternativa para que las líneas no se superpongan en el gráfico). Si agregué atan2 a esa gráfica, entonces sería la línea recta y = x. Entonces, me parece que podría haber varias respuestas, dependiendo de dónde se quiera poner la discontinuidad. Si uno realmente quiere replicar atan2 , la respuesta (en este género) sería

# Input: dx, dy: coordinates of a (difference) vector. # Output: a number from the range [-2 .. 2] which is monotonic # in the angle this vector makes against the x axis. # and with the same discontinuity as atan2 def pseudoangle(dx, dy): p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x if dy < 0: return p - 1 # -2 .. 0 increasing with x else: return 1 - p # 0 .. 2 decreasing with x

Esto significa que si el lenguaje que está utilizando tiene una función de signo, podría evitar la bifurcación devolviendo el signo (dy) (1-p), que tiene el efecto de poner una respuesta de 0 en la discontinuidad entre devolver -2 y +2. Y el mismo truco funcionaría con la metodología original de @ MvG, uno podría devolver el signo (dx) (p-1).

Actualización En un comentario a continuación, @MvG sugiere una implementación de C de una línea, a saber:

pseudoangle = copysign(1. - dx/(fabs(dx)+fabs(dy)),dy)

@MvG dice que funciona bien y me parece bien :-).

Muchos algoritmos (por ejemplo, la exploración de Graham ) requieren que los puntos o vectores se ordenen por su ángulo (tal vez como se ve desde algún otro punto, es decir, utilizando vectores de diferencia). Este orden es intrínsecamente cíclico, y cuando este ciclo se rompe para calcular valores lineales a menudo no importa mucho. Pero el valor del ángulo real tampoco importa mucho, siempre que se mantenga el orden cíclico. Así que hacer una llamada atan2 para cada punto podría ser un desperdicio. ¿Qué métodos más rápidos existen para calcular un valor que sea estrictamente monotónico en el ángulo, como atan2 es atan2 ? Dichas funciones aparentemente han sido llamadas "pseudoangle" por algunos.

Conozco una posible función de este tipo, que describiré aquí.

# Input: dx, dy: coordinates of a (difference) vector. # Output: a number from the range [-1 .. 3] (or [0 .. 4] with the comment enabled) # which is monotonic in the angle this vector makes against the x axis. def pseudoangle(dx, dy): ax = abs(dx) ay = abs(dy) p = dy/(ax+ay) if dx < 0: p = 2 - p # elif dy < 0: p = 4 + p return p

Entonces, ¿por qué funciona esto? Una cosa a tener en cuenta es que escalar todas las longitudes de entrada no afectará a la salida. Entonces, la longitud del vector (dx, dy) es irrelevante, solo importa su dirección. Concentrándonos en el primer cuadrante, por el momento podemos suponer que dx == 1 . Luego dy/(1+dy) crece monótonamente desde cero para dy == 0 a uno para infinito dy (es decir, para dx == 0 ). Ahora los otros cuadrantes tienen que ser manejados también. Si dy es negativo, entonces también lo es la p inicial. Así que para dx positivo ya tenemos un rango -1 <= p <= 1 monotónico en el ángulo. Para dx < 0 cambiamos el signo y agregamos dos. Eso da un rango de 1 <= p <= 3 para dx < 0 , y un rango de -1 <= p <= 3 en general. Si los números negativos son indeseables por alguna razón, se puede incluir la línea de comentarios de elif , que cambiará el cuarto cuadrante de -1…0 a 3…4 .

No sé si la función anterior tiene un nombre establecido y quién podría haberlo publicado primero. Lo obtuve hace bastante tiempo y lo copié de un proyecto a otro. Sin embargo, he encontrado occurrences de esto en la web, por lo que considero que este recorte es lo suficientemente público para su reutilización.

Hay una manera de obtener el rango [0 ... 4] (para ángulos reales [0 ... 2π]) sin introducir una distinción de caso adicional:

# Input: dx, dy: coordinates of a (difference) vector. # Output: a number from the range [0 .. 4] which is monotonic # in the angle this vector makes against the x axis. def pseudoangle(dx, dy): p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x if dy < 0: return 3 + p # 2 .. 4 increasing with x else: return 1 - p # 0 .. 2 decreasing with x

Lo más simple que se me ocurrió es hacer copias normalizadas de los puntos y dividir el círculo alrededor de ellos por la mitad a lo largo del eje x o y. Luego use el eje opuesto como un valor lineal entre el principio y el final del búfer superior o inferior (un búfer tendrá que estar en orden lineal inverso al ponerlo). Luego, puede leer el primer búfer y luego el segundo del búfer de manera lineal. sea ​​en el sentido de las agujas del reloj o segundo y primero en retroceso para el contrario.

Puede que esa no sea una buena explicación, así que puse un código en GitHub que usa este método para ordenar puntos con un valor de epsilion para dimensionar las matrices.

https://github.com/Phobos001/SpatialSort2D

Esto podría no ser bueno para su caso de uso porque está diseñado para el rendimiento en la representación de efectos gráficos, pero es rápido y simple (O (N) Complejidad). Si trabaja con cambios realmente pequeños en puntos o conjuntos de datos muy grandes (cientos de miles), esto no funcionará porque el uso de la memoria podría superar los beneficios del rendimiento.

Me gusta un poco la trigonometría, así que sé que la mejor forma de mapear un ángulo a algunos valores que usualmente tenemos es una tangente. Por supuesto, si queremos un número finito para no tener la molestia de comparar {signo (x), y / x}, se vuelve un poco más confuso.

Pero hay una función que mapea [1, + inf [to [1,0 [conocido como inverso], que nos permitirá tener un rango finito al cual mapearemos los ángulos. La inversa de la tangente es la cotangente bien conocida, por lo tanto, x / y (sí, es tan simple como eso).

Una pequeña ilustración, que muestra los valores de tangente y cotangente en un círculo unitario:

Verás que los valores son los mismos cuando | x | = | y |, y también ve que si coloreamos las partes que generan un valor entre [-1,1] en ambos círculos, podemos colorear un círculo completo. Para que este mapeo de valores sea continuo y monótono, podemos hacer dos cosas:

  • usa el opuesto de la cotangente para tener la misma monotonía que la tangente
  • agregue 2 a -cotan, para que los valores coincidan donde tan = 1
  • agregue 4 a la mitad del círculo (por ejemplo, debajo de x = -y diagonal) para que los valores se ajusten a la de las discontinuidades.

Eso da la siguiente función por partes, que es una función continua y monótona de los ángulos, con una sola discontinuidad (que es la mínima):

double pseudoangle(double dx, double dy) { // 1 for above, 0 for below the diagonal/anti-diagonal int diag = dx > dy; int adiag = dx > -dy; double r = !adiag ? 4 : 0; if (dy == 0) return r; if (diag ^ adiag) r += 2 - dx / dy; else r += dy / dx; return r; }

Tenga en cuenta que esto está muy cerca de los steve.hollasch.net/cgindex/math/fowler.html , con las mismas propiedades. Formalmente, pseudoangle(dx,dy) + 1 % 8 == Fowler(dx,dy)

Para hablar de rendimiento, es mucho menos ramificado que el código de Fowler (y, en general, una imo menos complicada). Compilada con -O3 en gcc 6.1.1, la función anterior genera un código de ensamblaje con 4 ramas, donde dos de ellas provienen de dy == 0 (una comprobando si ambos operandos están "desordenados", por lo tanto si dy era NaN , y el otro comprobando si son iguales).

Yo diría que esta versión es más precisa que otras, ya que solo utiliza las operaciones de preservación de la mantisa, hasta que se cambia el resultado al intervalo correcto. Esto debería ser especialmente visible cuando | x | << | y | o | y | >> | x |, entonces la operación | x | + | y ​​| Pierde bastante precisión.

Como puede ver en la gráfica, la relación ángulo-pseudoangulo también es muy cercana a la lineal.

Mirando de dónde vienen las ramas, podemos hacer las siguientes observaciones:

  • Mi código no se basa en los abs ni en la copysign , lo que lo hace parecer más autónomo. Sin embargo, jugar con bits de signo en valores de punto flotante es en realidad más bien trivial, ya que solo se voltea un bit separado (¡sin rama!), Así que esto es más una desventaja.

  • Además, otras soluciones propuestas aquí no comprueban si abs(dx) + abs(dy) == 0 antes de dividir por él, pero esta versión fallaría tan pronto como solo un componente (dy) sea 0, de modo que se lance en una rama (o 2 en mi caso).

Si elegimos obtener aproximadamente el mismo resultado (hasta errores de redondeo) pero sin sucursales, podríamos abusar de la firma y escribir:

double pseudoangle(double dx, double dy) { double s = dx + dy; double d = dx - dy; double r = 2 * (1.0 - copysign(1.0, s)); double xor_sign = copysign(1.0, d) * copysign(1.0, s); r += (1.0 - xor_sign); r += (s - xor_sign * d) / (d + xor_sign * s); return r; }

Pueden ocurrir errores más grandes que con la implementación anterior, debido a la cancelación en d o s si dx y dy están cerca en valor absoluto. No hay verificación de que la división por cero sea comparable con las otras implementaciones presentadas, y porque esto solo ocurre cuando dx y dy son 0.

Sólo tiene que utilizar una función de producto cruzado. La dirección en la que gire un segmento en relación con el otro dará un número positivo o negativo. No hay funciones trigonométricas ni división. Rápido y sencillo. Sólo buscalo en Google.

Si puede alimentar los vectores originales en lugar de los ángulos en una función de comparación al ordenar, puede hacer que funcione con:

  • Sólo una sola rama.
  • Solo comparaciones de punto flotante y multiplicaciones.

Evitar sumas y restas lo hace numéricamente mucho más robusto. En realidad, un doble siempre puede representar exactamente el producto de dos flotadores, pero no necesariamente su suma. Esto significa que para una entrada de precisión única puede garantizar un resultado perfecto y perfecto con poco esfuerzo.

Esta es básicamente share repetida para ambos vectores, con las ramas eliminadas y las divisiones multiplicadas. Devuelve un entero, con el signo que coincide con el resultado de la comparación (positivo, negativo o cero):

signed int compare(double x1, double y1, double x2, double y2) { unsigned int d1 = x1 > y1; unsigned int d2 = x2 > y2; unsigned int a1 = x1 > -y1; unsigned int a2 = x2 > -y2; // Quotients of both angles. unsigned int qa = d1 * 2 + a1; unsigned int qb = d2 * 2 + a2; if(qa != qb) return((0x6c >> qa * 2 & 6) - (0x6c >> qb * 2 & 6)); d1 ^= a1; double p = x1 * y2; double q = x2 * y1; // Numerator of each remainder, multiplied by denominator of the other. double na = q * (1 - d1) - p * d1; double nb = p * (1 - d1) - q * d1; // Return signum(na - nb) return((na > nb) - (na < nb)); }

agradable ... aquí hay una variante que devuelve -Pi, Pi como muchas funciones arctan2.

Nota de edición: cambié mi pseudocódigo a Python apropiado. Se cambió el orden de arg para que sea compatible con el módulo matemático Pythons atan2 (). Edit2 molesta más código para detectar el caso dx = 0.

def pseudoangle( dy , dx ): """ returns approximation to math.atan2(dy,dx)*2/pi""" if dx == 0 : s = cmp(dy,0) else:: s = cmp(dx*dy,0) # cmp == "sign" in many other languages. if s == 0 : return 0 # doesnt hurt performance much.but can omit if 0,0 never happens p = dy/(dx+s*dy) if dx < 0: return p-2*s return p

En esta forma, el error máximo es solo ~ 0.07 radianes para todos los ángulos. (por supuesto, omita el Pi / 2 si no le importa la magnitud).

Ahora, para la mala noticia, en mi sistema, usar python math.atan2 es aproximadamente un 25% más rápido. Obviamente, reemplazar un código interpretado simple no es mejor que una función compilada.