Ubicación de mayor densidad en una esfera

Tengo muchos puntos en la superficie de la esfera. ¿Cómo puedo calcular el área / punto de la esfera que tiene la mayor densidad de puntos? Necesito que esto se haga muy rápido. Si esto fuera un cuadrado, por ejemplo, supongo que podría crear una cuadrícula y luego dejar que los puntos voten qué parte de la cuadrícula es la mejor. He intentado transformar los puntos en coordenadas esféricas y luego hacer una cuadrícula, esto no funcionó bien, ya que los puntos alrededor del polo norte están cerca de la esfera pero distantes después de la transformación.

Gracias

Para agregar algunos otros esquemas alternativos a la mezcla: es posible definir un número de cuadrículas (casi) regulares en geometrías esféricas al refinar un poliedro inscrito.

La primera opción se llama una grilla icosaédrica, que es una triangulación de la superficie esférica. Al unir los centros de los triangularjs alrededor de cada vértice, también puede crear una grilla doble hexagonal basada en la triangulación subyacente:

rejilla icosaédrica

Otra opción, si no le gustan los triangularjs (y / o hexágonos) es la cuadrícula de esferas en cubos, formada por la subdivisión de las caras de un cubo inscrito y la proyección del resultado en la superficie esférica:

enter image description here

En cualquier caso, el punto importante es que las cuadrículas resultantes son casi regulares, por lo que para evaluar la región de mayor densidad en la esfera, simplemente puede realizar un análisis de estilo de histogtwig, contando el número de muestras por celda de cuadrícula.

Como han señalado varios comentaristas, para tener en cuenta la ligera irregularidad en la cuadrícula, es posible normalizar los recuentos del histogtwig al dividir por el área de cada celda de la cuadrícula. La densidad resultante se da luego como una medida “por unidad de área”. Para calcular el área de cada celda de la grilla hay dos opciones: (i) se puede calcular el área “plana” de cada celda, suponiendo que los bordes son líneas rectas, una aproximación que es probablemente bastante buena cuando la grilla es lo suficientemente densa, o (ii) puede calcular las áreas superficiales “verdaderas” mediante la evaluación de las integrales superficiales necesarias.

Si está interesado en realizar las consultas requeridas de “punto en celda” de manera eficiente, un enfoque es construir la cuadrícula como un árbol cuadrúpedo , comenzando con un poliedro inscrito y refinando sus caras en un árbol de caras secundarias. Para ubicar la celda adjunta, simplemente puede atravesar el árbol desde la raíz, que normalmente es una operación O(log(n)) .

Puede obtener información adicional sobre estos tipos de grillas aquí .

De hecho, no hay una razón real para dividir la esfera en una malla regular que no se superpone, intente esto:

  • particiona tu esfera en círculos semi-superpuestos

    mira aquí para generar puntos uniformemente distribuidos (tus centros circulares)

    Dispersión de n puntos uniformemente en una esfera

  • puede identificar los puntos en cada círculo muy rápido mediante un producto de puntos simple … realmente no importa si se cuentan dos puntos, el círculo con la mayor cantidad de puntos aún representa la densidad más alta

enter image description here

la implementación de mathematica

esto toma 12 segundos para analizar 5000 puntos. (y tomó aproximadamente 10 minutos para escribir)

  testcircles = { RandomReal[ {0, 1}, {3}] // Normalize}; Do[While[ (test = RandomReal[ {-1, 1}, {3}] // Normalize ; Select[testcircles , #.test > .9 & , 1] ) == {} ]; AppendTo[testcircles, test];, {2000}]; vmax = testcircles[[First@ Ordering[-Table[ Count[ (testcircles[[i]].#) & /@ points , x_ /; x > .98 ] , {i, Length[testcircles]}], 1]]]; 

enter image description here

Particione la esfera en regiones de área igual (delimitadas por paralelos y meridianos) como se describe en mi respuesta allí y cuente los puntos en cada región.

La relación de aspecto de las regiones no será uniforme (las regiones ecuatoriales serán más “cuadradas” cuando N~M , mientras que las regiones polares serán más alargadas). Esto no es un problema porque los diámetros de las regiones van a 0 a medida que N y M aumentan. La simplicidad computacional de este método prevalece sobre la mejor uniformidad de los dominios en las otras excelentes respuestas que contienen bellas imágenes.

Una modificación simple sería agregar dos regiones de “casquete polar” a las regiones N*M descritas en la respuesta vinculada para mejorar la estabilidad numérica (cuando el punto está muy cerca de un polo, su longitud no está bien definida). De esta forma, la relación de aspecto de las regiones está limitada.

Tratar puntos en una esfera como puntos 3D podría no ser tan malo.

Pruebe cualquiera de los siguientes:

  1. Seleccione k, realice una búsqueda aproximada de k-NN en 3D para cada punto en los datos o punto de interés seleccionado, luego pondere el resultado por su distancia hasta el punto de consulta. La complejidad puede variar para diferentes algoritmos k-NN aproximados.
  2. Cree una estructura de datos de partición de espacio como kd Tree, luego haga una consulta de conteo de rango aproximado (o exacta) con un scope de bola centrado en cada punto de los datos o punto de interés seleccionado. La complejidad es O (log (n) + epsilon ^ (- 3)) o O (epsilon ^ (- 3) * log (n)) para cada consulta de rango aproximado con algoritmos de última generación, donde épsilon es el umbral de error de rango wrt el tamaño de la bola que consulta. Para la consulta de rango exacto, la complejidad es O (n ^ (2/3)) para cada consulta.

Puede usar la proyección de Peters , que conserva las áreas.

Esto le permitirá contar eficientemente los puntos en una cuadrícula, pero también en una ventana deslizante (ventana Parzen de cuadro) utilizando el truco de la imagen integral .

Si entiendo correctamente, estás tratando de encontrar el punto denso en la esfera.

si los puntos son más densos en algún punto

  • Considera las coordenadas cartesianas y encuentra la media X, Y, Z de puntos

  • Encuentre el punto más cercano para significar X, Y, Z que está en la esfera (puede considerar el uso de coordenadas esféricas, solo extienda el radio al radio original).

Restricciones

  • Si la distancia entre la media X, Y, Z y el centro es menor que r / 2, entonces este algoritmo puede no funcionar como se desee.

No soy un maestro de las matemáticas, pero puede ser que pueda resolverlo analíticamente de la siguiente manera:

1. Cortar la coordenada

2.R = (Σ (n = 0. N = max) (Σ (m = 0. M = n) (1 / A ^ diff_in_consecative)) * ángulo) / Σangle

A = puede cualquier constante

Esto es solo un inverso de esta respuesta mía

simplemente invierta las ecuaciones de los vértices de la superficie de la esfera equidistante en el índice de la celda de la superficie. Ni siquiera intentes visualizar la celda diferente, luego circula o enloqueces. Pero si alguien realmente lo hace, por favor, publique el resultado aquí (y déjeme ahora)

Ahora solo crea un mapa de celda 2D y haz el cálculo de densidad en O(N) (como se hacen los histogtwigs) similar a lo que Darren Engwirda propone en su respuesta

Así es como se ve el código en C ++

 //--------------------------------------------------------------------------- const int na=16; // sphere slices int nb[na]; // cells per slice const int na2=na<<1; int map[na][na2]; // surface cells const double da=M_PI/double(na-1); // latitude angle step double db[na]; // longitude angle step per slice // sherical -> orthonormal void abr2xyz(double &x,double &y,double &z,double a,double b,double R) { double r; r=R*cos(a); z=R*sin(a); y=r*sin(b); x=r*cos(b); } // sherical -> surface cell void ab2ij(int &i,int &j,double a,double b) { i=double(((a+(0.5*M_PI))/da)+0.5); if (i>=na) i=na-1; if (i< 0) i=0; j=double(( b /db[i])+0.5); while (j< 0) j+=nb[i]; while (j>=nb[i]) j-=nb[i]; } // sherical <- surface cell void ij2ab(double &a,double &b,int i,int j) { if (i>=na) i=na-1; if (i< 0) i=0; a=-(0.5*M_PI)+(double(i)*da); b= double(j)*db[i]; } // init variables and clear map void ij_init() { int i,j; double a; for (a=-0.5*M_PI,i=0;i1.0) c=1.0; glColor3f(0.2,0.2,0.2); glCircle3D(x,y,z,x,y,z,0.45*da,0); // outline glColor3f(0.1,0.1,c ); glCircle3D(x,y,z,x,y,z,0.45*da,1); // filled by bluish color the more dense the cell the more bright it is } } //--------------------------------------------------------------------------- 

El resultado es así:

densidad superficial

así que ahora solo ve lo que está en la matriz del map[][] , puedes encontrar el mínimo / máximo de densidad global / local o lo que necesites … No olvides que el tamaño es map[na][nb[i]] donde i es el primer índice en conjunto. El tamaño de la cuadrícula está controlado por una constante y cm es solo la densidad de la escala de color …

[edit1] obtuvo la cuadrícula Quad, que es una representación mucho más precisa del mapeo usado

densidad superficial

esto es con na=16 los peores errores de redondeo están en los polos. Si quiere ser preciso, puede calcular la densidad de peso según el tamaño de la superficie de la celda. Para todas las células no polares, es simple quad. Para postes su triángulo ventilador (polígono regular)

Este es el código de cuadrícula:

 // draw cell quad grid (color is function of density) int i,j,ii,jj; double x,y,z,a,b,c,cm=1.0/10.0,mm=0.49,r=1.0; double dx=mm*da,dy; for (i=1;i1.0) c=1.0; glColor3f(0.2,0.2,0.2); glBegin(GL_LINE_LOOP); abr2xyz(x,y,z,a-dx,b-dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a-dx,b+dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a+dx,b+dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a+dx,b-dy,r); glVertex3d(x,y,z); glEnd(); glColor3f(0.1,0.1,c ); glBegin(GL_QUADS); abr2xyz(x,y,z,a-dx,b-dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a-dx,b+dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a+dx,b+dy,r); glVertex3d(x,y,z); abr2xyz(x,y,z,a+dx,b-dy,r); glVertex3d(x,y,z); glEnd(); } i=0; j=0; ii=i+1; dy=mm*db[ii]; ij2ab(a,b,i,j); c=map[i][j]; c=0.1+(c*cm); if (c>1.0) c=1.0; glColor3f(0.2,0.2,0.2); glBegin(GL_LINE_LOOP); for (j=0;j1.0) c=1.0; glColor3f(0.2,0.2,0.2); glBegin(GL_LINE_LOOP); for (j=0;j 

el mm es el tamaño de celda de la cuadrícula mm=0.5 es el tamaño de celda completa, menos crea un espacio entre las celdas

Considera usar un método geográfico para resolver esto. Las herramientas GIS, los tipos de datos geograficos en SQL, etc. manejan la curvatura de un esferoide. Puede que tenga que encontrar un sistema de coordenadas que use una esfera pura en lugar de un esferoide similar a la tierra si no está realmente modelando algo en la Tierra.

Para la velocidad, si tiene un gran número de puntos y desea la ubicación más densa de ellos, una solución tipo mapa de calor ttwig podría funcionar bien. Puede crear rásteres de baja resolución, luego hacer zoom en áreas de alta densidad y crear celdas de mayor resolución que le interesan.