I think to have found a little error in the file blob.cpp
responsible of one of the bugs reported under (job000189)
at http://www.povray.org/download/3.5-bugs.php
(exactly http://news.povray.org/3B991C90.D715B7A3@free.fr)

Look at the following snippet:

<CODE src=blob.cpp function=All_Blob_Intersections lines=~550ff>
    w = intervals[i+1].bound - l;

    newcoeffs[0] = coeffs[0] * w * w * w * w;
    newcoeffs[1] = (coeffs[1] + 4.0 * coeffs[0] * l) * w * w * w;
    newcoeffs[2] = (3.0 * l * (2.0 * coeffs[0] * l + coeffs[1]) + coeffs[2]) * w * w;
    newcoeffs[3] = (2.0 * l * (2.0 * l * (coeffs[0] * l + 0.75 * coeffs[1]) +
coeffs[2]) + coeffs[3]) * w;
    newcoeffs[4] = l * (l * (l * (coeffs[0] * l + coeffs[1]) + coeffs[2]) + coeffs[3])
+ coeffs[4];

    /* Calculate coefficients of corresponding bezier curve. [DB 10/94] */
</CODE>

I calculate those coefficients equating the two forms of a 4th order 
polynomial:

a0 t^4 + a1 t^3 + a2 t^2 + a3 t + a4 
=
   t^0 (1 - t)^4 * d0 + 
 4 t   (1 - t)^3 * d1 +
 6 t^2 (1 - t)^2 * d2 + 
 4 t^3 (1 - t)   * d3 +
   t^4 (1 - t)^0 * d4


Expanding the powers and reordering the rhs I obtain:

    (   d0) +
t   (-4 d0 +  4 d1) +
t^2 ( 6 d0 - 12 d1 +  6 d2) +
t^3 (-4 d0 + 12 d1 - 12 d2 + 4 d3) +
t^4 (   d0 -  4 d1 +  6 d2 - 4 d3 + d4)


Equating polynomial coefficients lead to the system:

a4 =    d0
a3 = -4 d0 +  4 d1
a2 =  6 d0 - 12 d1 +  6 d2
a1 = -4 d0 + 12 d1 - 12 d2 + 4 d3
a0 =    d0 -  4 d1 +  6 d2 - 4 d3 + d4

That's easily inverted

  d0 =   a4
4 d1 = 4 a4 +   a3
6 d2 = 6 a4 + 3 a3 +   a2
4 d3 = 4 a4 + 3 a3 + 2 a2 + a1
  d4 =   a4 +   a3 +   a2 + a1 + a0


Now comparing the code following the previous snippet with this 
solution I see three differences, in fact:

    d0 =   a4
    dk[0] = newcoeffs[4];

  4 d1 = 4 a4 +   a3
    dk[1] = newcoeffs[4] + 0.25 * newcoeffs[3];

  6 d2 = 6 a4 + 3 a3 +   a2
    dk[2] = newcoeffs[4] + 0.50 * (newcoeffs[3] + newcoeffs[2] / 12.0);

      but 1 / 6 !=         0.50 *                            1 / 12
      (and between 1 / 6 and 1 / 24 I prefer 1 / 6)


  4 d3 = 4 a4 + 3 a3 + 2 a2 + a1
    dk[3] = newcoeffs[4] + 0.50 * (0.375 * newcoeffs[3] + newcoeffs[2] + 0.125 *
newcoeffs[1]);

      but 3 / 4 !=         0.50 *  3 / 8
      (and between 3 / 4 and 3 / 16 I prefer 3 / 4)

      but 1 / 4 !=         0.50 *                                        1 / 8
      (and between 1 / 4 and 1 / 16 I prefer 1 / 4)

    d4 =   a4 +   a3 +   a2 + a1 + a0
    dk[4] = newcoeffs[4] + newcoeffs[3] + newcoeffs[2] + newcoeffs[1] + newcoeffs[0];


Now it seems that correcting these three errors the problem is fixed.

If someone has time and inclination to check what I did, I'll appreciate
and, as always, hints, comments and counter examples are welcome.

HTH Massimo

Post a reply to this message