Calculates view space 3D sphere extents on the screen

sphere_screen_extents.h

struct vec3f {float x, y, z;};
struct vec4f {float x, y, z, w;};
struct mat44f {vec4f x, y, z, w;};

//============================================================================
// sphere_screen_extents
//============================================================================
// Calculates the exact screen extents xyzw=[left, bottom, right, top] in
// normalized screen coordinates [-1, 1] for a sphere in view space. For
// performance, the projection matrix (v2p) is assumed to be setup so that
// z.w=1 and w.w=0. The sphere is also assumed to be completely in front
// of the camera.
// This is an optimized implementation of paper "2D Polyhedral Bounds of a
// Clipped Perspective-Projected 3D Sphere": http://jcgt.org/published/0002/02/05/paper.pdf
vec4f sphere_screen_extents(const vec3f &pos_, float rad_, const mat44f &v2p_)
{
  // calculate horizontal extents
  assert(v2p_.z.w==1 && v2p_.w.w==0);
  vec4f res;
  float rad2=rad_*rad_, d=pos_.z*rad_;
  float hv=sqrt(pos_.x*pos_.x+pos_.z*pos_.z-rad2);
  float ha=pos_.x*hv, hb=pos_.x*rad_, hc=pos_.z*hv;
  res.x=(ha-d)*v2p_.x.x/(hc+hb); // left
  res.z=(ha+d)*v2p_.x.x/(hc-hb); // right

  // calculate vertical extents
  float vv=sqrt(pos_.y*pos_.y+pos_.z*pos_.z-rad2);
  float va=pos_.y*vv, vb=pos_.y*rad_, vc=pos_.z*vv;
  res.y=(va-d)*v2p_.y.y/(vc+vb); // bottom
  res.w=(va+d)*v2p_.y.y/(vc-vb); // top
  return res;
}
//----------------------------------------------------------------------------