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

Casablanca

It's not really possible to answer as it stands. If you have a dot at (.5,.5) there's really high density in the square (.49999,.49999)(.500001,.500001). Do you have a minimum area or some other constraint? – Mantra

Try assigning weights to points depending on latitude. – Cameleer

aioobe: I guess I could infer some area constraints or something like that. For the moment I only have X points distributed on the sphere. Can I find some statistical distribution of the points on the sphere? ML? I could also infer a area constraints. How would I then proceed? – Casablanca

I think you will find it hard to produce an equal area grid. In Postgres/Postgis there is a geography data type that allows you to store points on a globe, add a spatial index, and then run queries using Intersects or Contains, which you could run with boxes or buffered points (circles). If you put this in a loop, you could get a grid density. Sorry, if this is a slightly off the wall suggestion. – Alrick

Cover the sphere with a recursive subdivision of an icosahedron and use the triangles thus formed. – Asha

@Tmyklebu how fast can i check which tringle has the greatest number of points? – Casablanca

@Johan: Linear time. Project your point to each of the xy-plane, the xz-plane, and the yz-plane and snap it to a grid. For each grid cell in each plane, make a list of the triangles that, when projected, hit that grid cell. If your grid coarseness is about the same as the size of a triangle, the smallest of the three cells will contain a constant number of triangles; use it. – Asha

Tmyklebu's suggestion, like Darren Engwirda's answer, introduces a bias that you have to correct. The correction involves computing some areas of small spherical triangles, which can be done using the Gauss-Bonnet formula and some trigonometry or approximated by computing the Jacobian of the projection map. I recommend using an approach which avoids introducing this bias instead. – Willful

You are wrong again in several places. The bias is not small. The difference between spherical area and the area in space can be quite large. If you Google the area of a spherical triangle, the first hit is Wolfram Mathworld's page giving the area in terms of the angular defect, the Gauss-Bonnet formula, as I said. This requires you to calculate some spherical angles, which are not the same as the 3d angles, and I'd bet most people would err in the rather involved calculations. So again, I recommend avoiding the approaches which lead to this bias if you actually want to implement this. – Willful

@DouglasZare: The Gauss-Bonnet formula is unnecessary; it's just linear algebra. If a spherical triangle on the unit sphere with aspect ratio at most 2 has area epsilon, the Euclidean triangle has area epsilon + O(epsilon^2), with a very small constant that I can't be bothered to work out. Please stop recommending approaches that are prone to give entirely wrong answers to the poster's question. – Asha

@tmyklebu: If computing spherical areas is so easy, check that your statement about the bias being small is wrong. The area can be distorted by a factor of 1.584 using an icosahedral subdivision with larger relative errors. I don't know where you came up with the epsilon+O(epsilon^2) but there is no such bound on the area of the triangle in space from the area of the spherical triangle. A triangle in space can have area 1 with a spherical area arbitrarily close to 0. The approaches by sds and george, among others, are much simpler and less error-prone than trying to do spherical trigonometry. – Willful

@DouglasZare: Bounds on area and aspect ratio imply the asymptotic bound I gave. The approach by sds does not work unmodified. The approach by george will involve a nontrivial spatial partitioning data structure. And, as I mentioned, the spherical trigonometry is straightforward and not particularly necessary. – Asha

Just consider Cartesian coordinates and find the mean X,Y,Z of points, it will be inside the sphere but you may find the closest point on sphere surface, that is your densepoint. More detailed answer is below as suggested answers. – Batty

Answer 1

To add some other, alternative schemes to the mix: it's possible to define a number of (almost) regular grids on sphere-like geometries by refining an inscribed polyhedron.

The first option is called an icosahedral grid, which is a triangulation of the spherical surface. By joining the centres of the triangles about each vertex, you can also create a dual hexagonal grid based on the underlying triangulation:

Another option, if you dislike triangles (and/or hexagons) is the cubed-sphere grid, formed by subdividing the faces of an inscribed cube and projecting the result onto the spherical surface:

In either case, the important point is that the resulting grids are almost regular -- so to evaluate the region of highest density on the sphere you can simply perform a histogram-style analysis, counting the number of samples per grid cell.

As a number of commenters have pointed out, to account for the slight irregularity in the grid it's possible to normalise the histogram counts by dividing through by the area of each grid cell. The resulting density is then given as a "per unit area" measure. To calculate the area of each grid cell there are two options: (i) you could calculate the "flat" area of each cell, by assuming that the edges are straight lines -- such an approximation is probably pretty good when the grid is sufficiently dense, or (ii) you can calculate the "true" surface areas by evaluating the necessary surface integrals.

If you are interested in performing the requisite "point-in-cell" queries efficiently, one approach is to construct the grid as a quadtree -- starting with a coarse inscribed polyhedron and refining it's faces into a tree of sub-faces. To locate the enclosing cell you can simply traverse the tree from the root, which is typically an O(log(n)) operation.

You can get some additional information regarding these grid types here.

Roundsman

The "cube" option does not seem to preserve areas, hence it will introduce a bias that needs to be corrected. – Propylite

The pictures are pretty, but the formulas are missing. How do you determine that a specific vector belongs to the specific domain? – Ribbonwood

@sds: In the cube-projection case, you project your points onto the cube. In the geodesic-tessellation case, you can still project your points onto the cube and then do O(1) point-in-triangle tests as I mentioned in my comments. The only thing really missing from this is the normalisation issue pointed out by Yves Daoust. – Asha

Answer 2

There is in fact no real reason to partition the sphere into a regular non-overlapping mesh, try this:

• partition your sphere into semi-overlapping circles

  see here for generating uniformly distributed points (your circle centers)

  Dispersing n points uniformly on a sphere

• you can identify the points in each circle very fast by a simple dot product..it really doesn't matter if some points are double counted, the circle with the most points still represents the highest density

mathematica implementation

this takes 12 seconds to analyze 5000 points. (and took about 10 minutes to write )

 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]]];

Oliver

This answer has the key: it doesn't matter if points are double counted! – Tjader

You're consigning yourself to either quadratically many dot products or some sort of spatial partitioning data structure with this approach. – Asha

there are a bunch of interesting answers here. The best approach is going to depend on details, such as how many points, what area to average over, etc. – Oliver

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.

Ribbonwood

Having the diameters of the regions tend to zero isn't really a strong enough guarantee; you need them to be bounded close enough to zero at a modest grid size since too fine a grid will essentially just pick out some happenstance distribution of noise. Separate tiling of the polar caps is a good idea if you choose to go this route, but at that stage you're closer to Darren Engwirda's second answer. – Asha

Answer 4

Treating points on a sphere as 3D points might not be so bad.

Try either:

1. Select k, do approximate k-NN search in 3D for each point in the data or selected point of interest, then weight the result by their distance to the query point. Complexity may vary for different approximate k-NN algorithms.
2. Build a space-partitioning data structure like k-d Tree, then do approximate (or exact) range counting query with a ball range centered at each point in the data or selected point of interest. Complexity is O(log(n) + epsilon^(-3)) or O(epsilon^(-3)*log(n)) for each approximate range query with state of the art algorithms, where epsilon is the range error threshold w.r.t. the size of the querying ball. For exact range query, the complexity is O(n^(2/3)) for each query.

Spineless

Answer 5

You can use the Peters projection, which preserves the areas.

This will allow you to efficiently count the points in a grid, but also in a sliding window (box Parzen window) by using the integral image trick.

Propylite

Interesting. What would projected squares on the grid correspond to on the sphere? Would they look similar? – Casablanca

Imagine a grid formed by equidistant meridians, and parallels with larger and larger spacing as you get closer to the poles. – Propylite

Answer 6

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

Entry

Answer 7

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.

Batty

Answer 8

I am not master of mathematics but may be it can solve by analytical way as:

1.Short the coordinate

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

A=may any constant

Apical

Answer 9

If you want a radial region of the greatest density, this is the robust disk covering problem with k = 1 and dist(a, b) = great circle distance (a, b) (see https://en.wikipedia.org/wiki/Great-circle_distance)

https://www4.comp.polyu.edu.hk/~csbxiao/paper/2003%20and%20before/PDCS2003.pdf

Adrea

Answer 10

Consider using a geographic method to solve this. GIS tools, geography data types in SQL, etc. all handle curvature of a spheroid. You might have to find a coordinate system that uses a pure sphere instead of an earthlike spheroid if you are not actually modelling something on Earth.

For speed, if you have large numbers of points and want the densest location of them, a raster heatmap type solution might work well. You could create low resolution rasters, then zoom to areas of high density and create higher resolution only cells that you care about.

Tinge