Acceleration of Disc primitive

Original:
static int Intersect_Disc (RAY *Ray, DISC *disc, DBL *Depth)
 {
  DBL t, u, v, r2, len;
  VECTOR P, D;

  Increase_Counter(stats[Ray_Disc_Tests]);

  /* Transform the point into the discs space */
  MInvTransPoint(P, Ray->Initial, disc->Trans);
  MInvTransDirection(D, Ray->Direction, disc->Trans);

  VLength(len, D);
  VInverseScaleEq(D, len);

  if (fabs(D[Z]) > EPSILON)
   {
    t = -P[Z] / D[Z];
    if (t >= 0.0)
     {
      u = P[X] + t * D[X];
      v = P[Y] + t * D[Y];
      r2 = Sqr(u) + Sqr(v);
      if ((r2 >= disc->iradius2) && (r2<= disc->oradius2))
       {
        *Depth = t / len;
        if ((*Depth > Small_Tolerance) && (*Depth < Max_Distance))
         {
          Increase_Counter(stats[Ray_Disc_Tests_Succeeded]);
          return (TRUE);
         }
       }
     }
   }
  return (FALSE);
 }
 

After optimization:
static int Intersect_Disc (RAY *Ray, DISC *disc, DBL *Depth)
 {
  DBL u, v, r2;
  VECTOR P, D;
  Increase_Counter(stats[Ray_Disc_Tests]);

  /* Transform the point into the discs space */
  MInvTransPoint(P, Ray->Initial, disc->Trans);
  MInvTransDirection(D, Ray->Direction, disc->Trans);

  if (fabs(D[Z]) > EPSILON)
   {
    *Depth = -P[Z] / D[Z];
    if ( *Depth > Small_Tolerance )
     {
      u = P[X] + *Depth * D[X];
      v = P[Y] + *Depth * D[Y];
      r2 = Sqr(u) + Sqr(v);
      if ( ( r2 >= disc->iradius2 ) && ( r2 <= 1.0 ) )
       {
        if ( *Depth < Max_Distance )
         {
          Increase_Counter(stats[Ray_Disc_Tests_Succeeded]);
          return (TRUE);
         }
       }
     }
   }
  return (FALSE);
 }

In above optimization we consider that:
- *Depth do not change their value after we normalize direction vector,
- Outer radius is hidden scale.

Same solution I was find in ART ( Another Ray Tracer ).
Unfortunatly I lost all information of that tracer.


If we change approach we have:

static int Intersect_Disc( RAY *Ray, DISC *Disc, DBL *Depth )
 {
  DBL u, v, r2;
  DBL NormalDotP, NormalDotDirection;
  VECTOR P, IPoint;

  Increase_Counter(stats[Ray_Disc_Tests]);

  VDot( NormalDotDirection, Disc->normal, Ray->Direction );
  if( fabs( NormalDotDirection ) < EPSILON ) return (FALSE);

  VSub( P, Ray->Initial, Disc->center );
  VDot( NormalDotP, Disc->normal, P );
  *Depth = - NormalDotP / NormalDotDirection;

  if ( *Depth <= Small_Tolerance ) return (FALSE);

  switch( Disc->Dominant_Axis )
   {
    case( X ): IPoint[ Y ] = P[ Y ] + *Depth * Ray->Direction[ Y ];
               IPoint[ Z ] = P[ Z ] + *Depth * Ray->Direction[ Z ];
               u = IPoint[ Y ] * Disc->u_axis[ X ] + IPoint[ Z ] * Disc->u_axis[ Y ];
               v = IPoint[ Y ] * Disc->v_axis[ X ] + IPoint[ Z ] * Disc->v_axis[ Y ];
     break;
    case( Y ): IPoint[ X ] = P[ X ] + *Depth * Ray->Direction[ X ];
               IPoint[ Z ] = P[ Z ] + *Depth * Ray->Direction[ Z ];
               u = IPoint[ Z ] * Disc->u_axis[ X ] + IPoint[ X ] * Disc->u_axis[ Y ];
               v = IPoint[ Z ] * Disc->v_axis[ X ] + IPoint[ X ] * Disc->v_axis[ Y ];
     break;
    case( Z ): IPoint[ X ] = P[ X ] + *Depth * Ray->Direction[ X ];
               IPoint[ Y ] = P[ Y ] + *Depth * Ray->Direction[ Y ];
               u = IPoint[ X ] * Disc->u_axis[ X ] + IPoint[ Y ] * Disc->u_axis[ Y ];
               v = IPoint[ X ] * Disc->v_axis[ X ] + IPoint[ Y ] * Disc->v_axis[ Y ];
     break;
   }

  r2 = Sqr(u) + Sqr(v);
  if ( ( Disc->iradius2 <= r2 ) && ( r2 <= 1.0 ) )
   {
    if ( (*Depth < Max_Distance) )
     {
      Increase_Counter(stats[Ray_Disc_Tests_Succeeded]);
      return (TRUE);
    }
   }
  return (FALSE);
 }
Looks like triangle intersect.

Tests & Source ( 41.873b )