...

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 / len;
  Intersection[i++].t = SIDE_HIT;
}

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

if ((t2 > Cone_Tolerance) && (t1 < Max_Distance) && (z >= 0.0) && (z <=
1.0))
{
  Intersection[i].d   = t2 / len;

...

This is a snippet from cones.cpp (lines ~ 220-235). In the second 'if'
I think the second test should be .. && (t2 < Max_Distance). And the
same typo is at line 291 in the same file.

The question concerns the divisions t1 / len  and  t2 /len. I think
that you should test for the validity of the intersection after the
division because otherwise you could discard some objects based on
a 'transformed' distance rather the real one.

The following example show the problem, you have 2 objects of the same
size and at the same distance from the camera. With a value of 100000
for the variable Scale you see only one, while if you use a value of
10000 you'll see both. The simple way to solve it, is to do the division
before the test. An other one could be to avoid the normalization of the
direction of the ray, it appears to me useless.

HTH Massimo


#declare Scale=100000;

camera { location <30*Scale, 0, 0>  look_at <0, 0, 0> }
light_source { <15*Scale, 0, 0>, rgb <1, 1, 1> }

cylinder {
  <0, 0, 0>, <-1, 0, 0> Scale
  pigment { color rgb <1, 1, 0> }
  scale <1, 5, 10>
  translate <0,0,10*Scale>
}
cylinder {
  <0, 0, 0>, <-1, 0, 0> 10*Scale
  pigment { color rgb <1, 0, 0> }
  scale <.1, .5, 1>
  translate <0,0,-10*Scale>
}

Post a reply to this message

Massimo Valentini wrote:


> 
> This is a snippet from cones.cpp (lines ~ 220-235). In the second 'if'
> I think the second test should be .. && (t2 < Max_Distance). And the
> same typo is at line 291 in the same file.
> 


I seems to agree with you, at least on the t1 vs t2 issue.


> The question concerns the divisions t1 / len  and  t2 /len. I think
> that you should test for the validity of the intersection after the
> division because otherwise you could discard some objects based on
> a 'transformed' distance rather the real one.
> 


I do not see why there is a limitation at Max_Distance,
I would simply remove that part of the test (It looks like an old 
optimisation: if it is enough away, remove it... it's bad, IMHO).


[...]

> An other one could be to avoid the normalization of the
> direction of the ray, it appears to me useless.

It is NOT useless. When CSG operation/object (union & merge, at least) 
will ask its component, it will orders the intersections by depth order 
before answering to the upper level. Thus, putting back the depth from 
object space to global space is mandatory (otherwise, strange things 
will happens).

Same goes when more than a single cone/cylinder get intersected by the 
ray, the engine need to order them correctly, in the global space!

Post a reply to this message

"Le Forgeron" ha scritto
:
: > An other one could be to avoid the normalization of the
: > direction of the ray, it appears to me useless.
:
: It is NOT useless. When CSG operation/object (union & merge, at least)
: will ask its component, it will orders the intersections by depth order
: before answering to the upper level. Thus, putting back the depth from
: object space to global space is mandatory (otherwise, strange things
: will happens).
:
: Same goes when more than a single cone/cylinder get intersected by the
: ray, the engine need to order them correctly, in the global space!
:

I try to show what I mean. You are first scaling the direction

<CODE src=cones.cpp lines=197-200>
MInvTransDirection(D, Ray->Direction, Cone->Trans);

VLength(len, D);
VInverseScaleEq(D, len);
</CODE>

This last line is a macro for a vector/scalar division:
D = D / len

Doing so results in a scaling of the coefficients (a & b)
of the 2nd order equation you are going to solve:

<CODE src=cones.cpp lines=206-213>
a = D[X] * D[X] + D[Y] * D[Y];

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

  c = P[X] * P[X] + P[Y] * P[Y] - 1.0;
</CODE>


And now you are required to scale the solution:

<CODE src= cones.cpp lines=220-228>
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 / len;
</CODE>


Suppose not to scale D, I name this not scaled direction D`
and D the direction that's actually used. It is:

D` = D * len

The resulting equation would be

a` * t`^2 + b` * t` + c = 0

where

a` = D`[X] * D`[X] + D`[Y] * D`[Y]
   = (D[X] * D[X] + D[Y] * D[Y]) * len ^ 2
   = a * len ^ 2

and b` = b * len

and whose solutions are t1,2` = t1,2 / len.

a` * t`^2 + b` * t` + c =
a * len^2 * (t / len)^2 + b * len * (t / len) + c =
a * t ^ 2 + b * t + c = 0

Exactly the previous solution only already scaled, so ready to
be used in global space sorting, testing, whatever necessary.

And

z` = P[Z] + t` * D`[Z]
   = P[Z] + (t / len) * (D[Z] * len)
   = P[Z] + t * D[Z] = z

The same z ready to be tested against 0 and 1.

So avoiding the normalization of the direction only saves a sqrt
(no need to calculate len) and some multiplications and divisions
but gives you the same result this function is returning.

The same reasoning is valid for the cone intersections.

If you want to try what I'm suggesting, assign 1.0 to len where this
is declared, and then comment out the computation of len.

<CODE lines=188-200>
  DBL a, b, c, z, t1, t2, len = 1.0;//<------
  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);

//  VLength(len, D);// <------------------------
//  VInverseScaleEq(D, len); //<----------------


HTH Massimo

Post a reply to this message

Massimo Valentini wrote:

> "Le Forgeron" ha scritto
> :
> : > An other one could be to avoid the normalization of the
> : > direction of the ray, it appears to me useless.
> :
> : It is NOT useless. When CSG operation/object (union & merge, at least)
> : will ask its component, it will orders the intersections by depth order
> : before answering to the upper level. Thus, putting back the depth from
> : object space to global space is mandatory (otherwise, strange things
> : will happens).
> :
> : Same goes when more than a single cone/cylinder get intersected by the
> : ray, the engine need to order them correctly, in the global space!
> :
> 
> I try to show what I mean. You are first scaling the direction
> 


Ok... I understand better what you wanted.
I first thought that you wanted only to remove the division of the depth.

So, far, your suggestion looks ok, and it might deserves a study of the 
code, as well as experiments (just to check there is no adverse effect 
due to a precision problem...probably not, but I'm afraid of sturm 
solver effects). The study should also find out if there is any speed-up 
(or slow down) with the new code. (and by how much).

Post a reply to this message

"Massimo Valentini" ha scritto 
: I try to show what I mean. You are first scaling the direction
: 
: <CODE src=cones.cpp lines=197-200>
: MInvTransDirection(D, Ray->Direction, Cone->Trans);
: 
: VLength(len, D);
: VInverseScaleEq(D, len);
: </CODE>
: 
: This last line is a macro for a vector/scalar division:
: D = D / len
: 
: Doing so results in a scaling of the coefficients (a & b)
: of the 2nd order equation you are going to solve:
: 
: <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. 

So you may either remove the condition for intersctions over 
Max_Distance, or shift the divisions t1,2/len before the test against 
Max_Distance. 

Massimo

Post a reply to this message

Massimo Valentini wrote:

> "Massimo Valentini" ha scritto 

> 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. 
> 


The (a > EPSILON) is used to avoid numerical error (division by 0 and 
alike, division by very small value gives a great inaccuracy)


EPSILON testing is also used later to avoid reflecting on the same 
object forever
(First ray hit the object, then you compute a ray from the hit point...
if EPSILON is not used, due to numerical error, the new ray may find 
again the same point [well, mathemacally, it is the same, but 
numerically it may be a little different so that testing would not 
detect it... hence the EPSILON : when the eye is near enough the 
surface, the surface should be removed... Epsilon is an arbitrary value,
so using EPSILON or EPSILON*len^2 should not be of too much concern for 
reasonnable value of len {ok, I might be contradicted by a scene...})

Sidenote: (a > EPSILON) and (t1 > EPSILON) show no correlation on the 
value of EPSILON, so there is probably no problem at all!


> So you may either remove the condition for intersctions over 
> Max_Distance, or shift the divisions t1,2/len before the test against 
> Max_Distance. 
> 


Test against Max_distance is unrelated with Epsilon value, and is any 
way a candidate for removal...

Post a reply to this message

"Le Forgeron" ha scritto 
: Massimo Valentini wrote:
: 
: > 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. 
: > 
: 
: 
: The (a > EPSILON) is used to avoid numerical error (division by 0 and 
: alike, division by very small value gives a great inaccuracy)
: 

This is what I thought at first. But this condition has a geometrical 
meaning too. In fact 'a' is the square of the sin of the angle formed 
by the ray and the cylindrical axis, when calculated with ||D|| (length
of D, the norm) == 1, while when D isn't normalized it is the sin of the 
same angle multiplied by the square of len . 

Obviously, if the ray is (close to be) parallel to the axis, the solution 
is either impossible or indeterminate, the ray does not intersect the 
cylindrical surface or it lies entirely on it.

Given the fact that I'm not aware of the range of variability of len 
after transformation, it appears now a little more complex to evaluate
whether the removal of the normalization has an impact only on some 
border case (tiny cylinders in kilometers wide panoramas) or it's 
influencing some real and well looking scene.

I was thinking to suggest a simple fix to a border case strangeness 
with no impact on the great majority of the scenes, but now it appears 
that a little more testing is necessary for evaluating both the 
performance and the rendering impact. 

Massimo

Post a reply to this message

"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

Massimo Valentini wrote:

> "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.


NO MORE INLINE here please (even binary... ).

I made also some tests based on 3.1 sources (I'm not yet ready for 3.5) 
and get also a feeling of speed, but I forgot that I updated my compiler 
in the meantime, so the result are meaningless (until I revert to the 
old code and compile again...)

Mathematically, the correction is ok.
correcting the EPSILON factor is IMHO useless, just keep it "(a>EPSILON)".

Post a reply to this message

In article <3DA### [at] freefr>, Le Forgeron <jgr### [at] freefr> 
wrote:

> NO MORE INLINE here please (even binary... ).

Huh? It is attachments that are forbidden outside the binary groups, 
in-line code like this is fine. This was a bit long, but not 
unacceptable in my opinion.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: chr### [at] tagpovrayorg
http://tag.povray.org/

Post a reply to this message