Location of highest density on a sphere

Question

I have a lot of points on the surface of the sphere. How can I calculate the area/spot of the sphere that has the largest point density? I need this to be done very fast. If this was a square for example I guess I could create a grid and then let the points vote which part of the grid is the best. I have tried with transforming the points to spherical coordinates and then do a grid, both this did not work well since points around north pole are close on the sphere but distant after the transform.

Thanks

Answer 1

This is really just an inverse of this answer of mine

just invert the equations of equidistant sphere surface vertexes to surface cell index. Don't even try to visualize the cell different then circle or you go mad. But if someone actually do it then please post the result here (and let me now)

Now just create 2D cell map and do the density computation in O(N) (like histograms are done) similar to what Darren Engwirda propose in his answer

This is how the code looks like in 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;i<na;i++,a+=da)
        {
        nb[i]=ceil(2.0*M_PI*cos(a)/da);     // compute actual circle cell count
        if (nb[i]<=0) nb[i]=1;
        db[i]=2.0*M_PI/double(nb[i]);       // longitude angle step
        if ((i==0)||(i==na-1)) { nb[i]=1; db[i]=1.0; }
        for (j=0;j<na2;j++) map[i][j]=0;    // clear cell map
        }
    }
//---------------------------------------------------------------------------
// this just draws circle from point x0,y0,z0 with normal nx,ny,nz and radius r
// need some vector stuff of mine so i did not copy the body here (it is not important)
void glCircle3D(double x0,double y0,double z0,double nx,double ny,double nz,double r,bool _fill);
//---------------------------------------------------------------------------
void analyse()
    {
    // n is number of points and r is just visual radius of sphere for rendering
    int i,j,ii,jj,n=1000;
    double x,y,z,a,b,c,cm=1.0/10.0,r=1.0;
    // init
    ij_init();      // init variables and map[][]
    RandSeed=10;    // just to have the same random points generated every frame (do not need to store them)
    // generate draw and process some random surface points
    for (i=0;i<n;i++)
        {
        a=M_PI*(Random()-0.5);
        b=M_PI* Random()*2.0 ;
        ab2ij(ii,jj,a,b);       // cell corrds
        abr2xyz(x,y,z,a,b,r);   // 3D orthonormal coords
        map[ii][jj]++; // update cell density
        // this just draw the point (x,y,z) as line in OpenGL so you can ignore this
        double w=1.1;   // w-1.0 is rendered line size factor
        glBegin(GL_LINES);
        glColor3f(1.0,1.0,1.0); glVertex3d(x,y,z);
        glColor3f(0.0,0.0,0.0); glVertex3d(w*x,w*y,w*z);
        glEnd();
        }

    // draw cell grid (color is function of density)
    for (i=0;i<na;i++)
     for (j=0;j<nb[i];j++)
        {
        ij2ab(a,b,i,j); abr2xyz(x,y,z,a,b,r);
        c=map[i][j]; c=0.1+(c*cm); if (c>1.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
        }
    }
//---------------------------------------------------------------------------

The result looks like this:

so now just see what is in the map[][] array you can find the global/local min/max of density or whatever you need... Just do not forget that the size is map[na][nb[i]] where i is the first index in array. The grid size is controlled by na constant and cm is just density to color scale ...

[edit1] got the Quad grid which is far more accurate representation of used mapping

this is with na=16 the worst rounding errors are on poles. If you want to be precise then you can weight density by cell surface size. For all non pole cells it is simple quad. For poles its triangle fan (regular polygon)

This is the grid draw code:

// 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;i<na-1;i++)    // ignore poles
 for (j=0;j<nb[i];j++)
    {
    dy=mm*db[i];
    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);
    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;j<nb[ii];j++) { ij2ab(a,b,ii,j); abr2xyz(x,y,z,a-dx,b-dy,r); glVertex3d(x,y,z); }
glEnd();
glColor3f(0.1,0.1,c  );
glBegin(GL_TRIANGLE_FAN);                 abr2xyz(x,y,z,a   ,b   ,r); glVertex3d(x,y,z);
for (j=0;j<nb[ii];j++) { ij2ab(a,b,ii,j); abr2xyz(x,y,z,a-dx,b-dy,r); glVertex3d(x,y,z); }
glEnd();
i=na-1; 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;j<nb[ii];j++) { ij2ab(a,b,ii,j); abr2xyz(x,y,z,a-dx,b+dy,r); glVertex3d(x,y,z); }
glEnd();
glColor3f(0.1,0.1,c  );
glBegin(GL_TRIANGLE_FAN);                 abr2xyz(x,y,z,a   ,b   ,r); glVertex3d(x,y,z);
for (j=0;j<nb[ii];j++) { ij2ab(a,b,ii,j); abr2xyz(x,y,z,a-dx,b+dy,r); glVertex3d(x,y,z); }
glEnd();

the mm is the grid cell size mm=0.5 is full cell size , less creates a space between cells

Answer 2

If I understand correctly, you are trying to find the densepoint on sphere.

if points are denser at some point

• Consider Cartesian coordinates and find the mean X,Y,Z of points

• Find closest point to mean X,Y,Z that is on sphere (you may consider using spherical coordinates, just extend the radius to original radius).

Constraints

• If distance between mean X,Y,Z and the center is less than r/2, then this algorithm may not work as desired.

Answer 3

Partition the sphere into equal-area regions (bounded by parallels and meridians) as described in my answer there and count the points in each region.

The aspect ratio of the regions will not be uniform (the equatorial regions will be more "squarish" when N~M, while the polar regions will be more elongated). This is not a problem because the diameters of the regions go to 0 as N and M increase. The computational simplicity of this method trumps the better uniformity of domains in the other excellent answers which contain beautiful pictures.

One simple modification would be to add two "polar cap" regions to the N*M regions described in the linked answer to improve the numeric stability (when the point is very close to a pole, its longitude is not well defined). This way the aspect ratio of the regions is bounded.