Spherical geometry optimisations

I am currently working on recoding proper classes for managing spherical geometry primitives in Stellarium. One of the most used primitive is the so called spherical cap showed in the figure. A spherical cap (sometimes called half space) is a disk-like region of the sphere defined by a direction unit vector n and an aperture angle ⍺.

The most common operation that we want to perform on a spherical cap is to check if a point A of the sphere lies inside the cap. This is pretty easy to do by checking whether the dot product n·OA is smaller than cos(⍺).

The C++ code below implements this. Note that to avoid computing a CPU expensive cosine each time we perform the contains operation, we pre-compute cos(⍺) and store it as a class member d.

class SphericalCap
{
    bool contains(const Vec3d &v) const {return (v*n>=d);}
    Vec3d n;
    double d;
};

Another useful operation we want to perform often is to check whether two spherical caps with direction vectors n1 and n2 and aperture ⍺1 and ⍺2 intersect. The straightforward implementation is to compute the angle θ between n1 and n2 and check whether it is smaller than ⍺1 + ⍺2.

bool intersects(const SphericalCap& other) const
{
    const double theta = n.angle(other.n);
    return theta < acos(d) + acos(other.d);
}

Unfortunately, because we store only the cosine of ⍺1 and ⍺2, this solution involves 2 calls to the acos function, plus 1 call to the angle() function which itself involves 1 acos and 1 sqrt. Needless to say that it's not lightning fast..

After some thoughts I finally came up with a solution involving only basic arithmetics:

First, in the case ⍺1+⍺2 >= 180°, (i.e. if d1+d2<=0) the 2 caps necessarily intersect so we don't need to think further. Let's now treat the other cases:

The 2 caps intersect only and only if:

θ <= ⍺1+⍺2         (1)

Because both members of (1) are < 180° (the angle θ between n1 and n2 can obviously not be > 180°), we can transform (1) into:

cos (θ) >= cos(⍺1+⍺2)

which with basic trigo can be transformed into:

L <= sin ⍺1 sin ⍺2 , with L = cos ⍺1 cos ⍺2 - cos θ     (2)

In our case L is known and given by d1*d2-n1·n2. Because ⍺1 and ⍺2 are bounded between 0° and 180°, sin ⍺1 sin ⍺2 is bounded between 0 and 1. Therefore if L<=0, (2) is verified and the 2 caps intersect. Similarily, if L>1, (2) is not verified and the 2 caps don't intersect. In the last case where 0<L<=1 we can transform (2) into:

L² <= sin²⍺1 sin²⍺2

which with basic trigo formulae can be transformed into:

L² <= (1-cos²⍺1)(1-cos²⍺2)         (3)

In our case cos²⍺1 and cos²⍺2 are know and given by d1² and d2² so we can easily check that (3) is verified or not.

We thus have a solution to each cases, and none of them involve any CPU expensive operation!! In C++ code, this reasoning translates into:

bool intersects(const SphericalCap& other) const
{
    const double L = d*other.d - n*other.n;
    return d+other.d<=0. || L<=0. || (L<=1. && L*L <= (1.-d*d)*(1.-other.d*other.d));
}