"Massimo Valentini" ha scritto
: : 
: : <CODE src=cones.cpp lines=206-213>
: : a = D[X] * D[X] + D[Y] * D[Y];
: : 
: : if (a > EPSILON)
: : {
: 
: This solution is not equivalent to the original code.
: 
: This condition is different whether you scale the direction 
: or not. May be there could be a way to avoid the square root (e.g.
: expressing this condition if (a > EPSILON*len^2)) but it'd require
: some work analyzing the code for both cones and cylinders. 
: 

It seems to be worth the effort, I patched the function intersect_cone 
in order to avoid the normalization of the ray direction and with some 
simple tests, based on scenes constituted mostly by cylinders, this 
patch saves a tenth of the total time and the image is the same produced 
by the original povray. 

Comments, hints and counter-examples are appreciated.

Massimo


Here's the patch and an example I used to time the differences:

static int intersect_cone(RAY *Ray, CONE *Cone, CONE_INT *Intersection)
{
  int i = 0;
  DBL a, b, c, z, t1, t2, squared_len;
  DBL d;
  VECTOR P, D;

  Increase_Counter(stats[Ray_Cone_Tests]);

  /* Transform the ray into the cones space */

  MInvTransPoint(P, Ray->Initial, Cone->Trans);
  MInvTransDirection(D, Ray->Direction, Cone->Trans);

  a = D[X] * D[X] + D[Y] * D[Y];
  squared_len = a + D[Z] * D[Z];

  if (Test_Flag(Cone, CYLINDER_FLAG))
  {
    /* Solve intersections with a cylinder */

    if (a > EPSILON * squared_len)
    {
      b = P[X] * D[X] + P[Y] * D[Y];

      c = P[X] * P[X] + P[Y] * P[Y] - 1.0;

      d = b * b - a * c;

      if (d >= 0.0)
      {
        d = sqrt(d);

        t1 = (-b + d) / a;
        t2 = (-b - d) / a;

        z = P[Z] + t1 * D[Z];

        if ((t1 > Cone_Tolerance) && (t1 < Max_Distance) && (z >= 0.0) && (z <= 1.0))
        {
          Intersection[i].d   = t1;
          Intersection[i++].t = SIDE_HIT;
        }

        z = P[Z] + t2 * D[Z];

        if ((t2 > Cone_Tolerance) && (t2 < Max_Distance) && (z >= 0.0) && (z <= 1.0))
        {
          Intersection[i].d   = t2;
          Intersection[i++].t = SIDE_HIT;
        }
      }
    }
  }
  else
  {
    /* Solve intersections with a cone */

    a -= D[Z] * D[Z];

    b = D[X] * P[X] + D[Y] * P[Y] - D[Z] * P[Z];

    c = P[X] * P[X] + P[Y] * P[Y] - P[Z] * P[Z];

    if (fabs(a) < EPSILON * squared_len)
    {
      if (fabs(b) > EPSILON * sqrt(squared_len))
      {
        /* One intersection */

        t1 = -c / (2 * b);

        z = P[Z] + t1 * D[Z];

        if ((t1 > Cone_Tolerance) && (t1 < Max_Distance) && (z >= Cone->dist) && (z <=
1.0))
        {
          Intersection[i].d   = t1;
          Intersection[i++].t = SIDE_HIT;
        }
      }
    }
    else
    {
      /* Check hits against the side of the cone */

      d = b * b - a * c;

      if (d >= 0.0)
      {
        d = sqrt(d);

        t1 = (-b - d) / a;
        t2 = (-b + d) / a;

        z = P[Z] + t1 * D[Z];

        if ((t1 > Cone_Tolerance) && (t1 < Max_Distance) && (z >= Cone->dist) && (z <=
1.0))
        {
          Intersection[i].d   = t1;
          Intersection[i++].t = SIDE_HIT;
        }

        z = P[Z] + t2 * D[Z];

        if ((t2 > Cone_Tolerance) && (t2 < Max_Distance) && (z >= Cone->dist) && (z <=
1.0))
        {
          Intersection[i].d   = t2;
          Intersection[i++].t = SIDE_HIT;
        }
      }
    }
  }

  if (Test_Flag(Cone, CLOSED_FLAG) && (fabs(D[Z]) > EPSILON))
  {
    d = (1.0 - P[Z]) / D[Z];

    a = (P[X] + d * D[X]);

    b = (P[Y] + d * D[Y]);

    if (((Sqr(a) + Sqr(b)) <= 1.0) && (d > Cone_Tolerance) && (d < Max_Distance))
    {
      Intersection[i].d   = d;
      Intersection[i++].t = CAP_HIT;
    }

    d = (Cone->dist - P[Z]) / D[Z];

    a = (P[X] + d * D[X]);

    b = (P[Y] + d * D[Y]);

    if ((Sqr(a) + Sqr(b)) <= (Test_Flag(Cone, CYLINDER_FLAG) ? 1.0 : Sqr(Cone->dist))
      && (d > Cone_Tolerance) && (d < Max_Distance))
    {
      Intersection[i].d   = d;
      Intersection[i++].t = BASE_HIT;
    }
  }

  if (i)
  {
    Increase_Counter(stats[Ray_Cone_Tests_Succeeded]);
  }

  return (i);
}


/*************************************/
#macro Dr(F,D) #local N = vlength(D);#declare D = D / N;
#while (N >= 0) cylinder {z,-z,0.3 translate F+D*N}
#declare N=N-1;#end#end

difference {
  cylinder { z*.09, z*.01, 9 }
  Dr(<-6, -3, 0>, 6*y)
  Dr(<-6, 3, 0>, <1.6, -2, 0>)
  Dr(<-3, 3, 0>,<-1.6, -2, 0>)
  Dr(<-3, -3, 0>, 6*y)
  Dr(<-2, -3, 0>, 5*y)
  Dr(<-1.75, 3, 0>, 3*x)
  Dr(<2, -3, 0>, 5*y)
  Dr(<-2, 0, 0>, 4*x)
  Dr(<3, -3.2, 0>, 1.1*<3, 6, 0>)
  Dr(<6, -3.2, 0>, 1.1*<-3, 6, 0>)
  pigment { color rgb <1, 0, 0> }
}
cylinder { <0,-100, 110> <0,100,110> 100 pigment {rgb <0,0,1>} }
light_source { <0, 0, -10>, rgb <1, 1, 1> }
camera {location <12, -15, -1.75> look_at 10*z }

Post a reply to this message