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

}

## No comments:

## Post a Comment