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

Massimo Valentini wrote:

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


You mean that you understood that optimisation code ? Greats if true!!

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



I have some time on this weekend, so I will have a look at it.
(My current suggestion for correction is to remove the optimisation,
which is really inappropriate if you are correct).
More on monday...

Post a reply to this message

"Le Forgeron" ha scritto 
: 
: You mean that you understood that optimisation code ? Greats if true!!
: 

What I understood is this:

You've got a polynomial and are interested in roots that lie
in the interval (l l+w):

b0 x^4 + b1 x^3 + b2 x^2 + b3 x + b4 = 0

First you change the coordinate with a transformation of this kind:

  x = l + t w    <--> t = (x - l) / w

so the interesting interval becomes (0 1), the substitution is trivial
and in the code it is correct AFAICS.

So you obtain a polynomial in t

a0 t^4 + a1 t^3 + a2 t^2 + a3 t + a4 = 0

you could already say that if all ajs are non negative (non positive)
then there is no root in (0 1), it's the sum of non negative () 
quantities..., but this is not really meaningful, since happens often
that one aj has a different sign from the others, but there is no root 
in that interval.

So he transforms the polynomial to the Bezier form, and the same 
reasoning holds (sum of non negative). And expressed in this form 
it is meaningful since when skipping the search the execution times 
decrease significantly.

I did not consider the case of all ajs == 0 or what happens in zero, 
one. But I'm confident in who wrote originally the code.

Massimo

Post a reply to this message

Le Forgeron wrote:
> 
> Massimo Valentini wrote:

> 
> > If someone has time and inclination to check what I did, I'll appreciate
> > and, as always, hints, comments and counter examples are welcome.
> 
> I have some time on this weekend, so I will have a look at it.
> (My current suggestion for correction is to remove the optimisation,
> which is really inappropriate if you are correct).
> More on monday...


No need to wait that long.
It's ok with all the bogus blob scene I had.

With the correction, it's faster for rendering than removing the optimisation, and 
there is far less need for sturm (which allow to keep more speed too!!!).

I would suggest that the correction be incorporated in the official code.

-- 
Non Sine Numine
http://grimbert.cjb.net/
Etiquette is for those with no breeding;
fashion for those with no taste.

Post a reply to this message

> I think to have found a little error in the file blob.cpp


Wow, this is a complicated post. While I may be able to understand it if I
slowly read all the way through it, would you care to give a simplified
explanation of what was wrong in the code? Was it just a typo? Was something
assumed to work that didn't? Just out of curiosity.

And was this the only bug in the blob code? That'd be nice if it were =)

 - Slime
[ http://www.slimeland.com/ ]

Post a reply to this message

"Slime" ha scritto 
: > I think to have found a little error in the file blob.cpp
: 
: 
: Wow, this is a complicated post. While I may be able to understand it if I
: slowly read all the way through it, would you care to give a simplified
: explanation of what was wrong in the code? Was it just a typo? Was something
: assumed to work that didn't? Just out of curiosity.

You may call it a typo. I think someone trying to simplify a computation
(division -> multiplication of the reciprocal or similar) made a double
division instead of a division and a multiplication, so three numbers 
ended up to be a fourth of their correct value.

The lines to modify are around 545 in the file blob.cpp (if you need more 
context, look at my previous post or the bug report). 
Substitute for example:

dk[2] = newcoeffs[4] + 0.50 * (newcoeffs[3] + newcoeffs[2] / 12.0);
dk[3] = newcoeffs[4] + 0.50 * (0.375 * newcoeffs[3] + newcoeffs[2] + 0.125 *
newcoeffs[1]);

with:

dk[2] = newcoeffs[4] + 0.50 * (newcoeffs[3] + newcoeffs[2] / 3.0);
dk[3] = newcoeffs[4] + 0.50 * (1.5 * newcoeffs[3] + newcoeffs[2] + 0.5 *
newcoeffs[1]);


: And was this the only bug in the blob code? That'd be nice if it were =)
: 

It is the only one I found. But it does not affect other scenes 
reported as bugs.

Massimo

Post a reply to this message

Massimo Valentini wrote:

> "Slime" ha scritto 

> 
> : And was this the only bug in the blob code? That'd be nice if it were =)
> : 
> 
> It is the only one I found. But it does not affect other scenes 
> reported as bugs.


Well, it might not correct the Media & Blob (I did not test that), but 
at least, it perfectly corrects the cuts that have been present in blob 
since a long time (and that started job 189 during the release-beta 
periode).
Moreover, your correction is mathematically explained and make sense.

I'm just wondering if any people of the Pov-team reads this newsgroup 
(and will take the solution from here) or if you need to make a formal 
posting somewhere else (p.bugreport is moderated, but now that you have 
discussed here, the moderator should agree to a post... otherwise, I'm 
afraid of a post in p.general, just to get the attention of the Pov-team 
members...)

Post a reply to this message

Le Forgeron wrote:

> I'm just wondering if any people of the Pov-team reads this newsgroup
> (and will take the solution from here) or if you need to make a formal
> posting somewhere else (p.bugreport is moderated, but now that you have
> discussed here, the moderator should agree to a post... otherwise, I'm
> afraid of a post in p.general, just to get the attention of the Pov-team
> members...)

I will check into this.

-- 
Ken Tyler

Post a reply to this message

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

> I'm just wondering if any people of the Pov-team reads this newsgroup
> (and will take the solution from here) or if you need to make a formal
> posting somewhere else (p.bugreport is moderated, but now that you have
> discussed here, the moderator should agree to a post... otherwise, I'm
> afraid of a post in p.general, just to get the attention of the Pov-team
> members...)

No patch has been posted.  Just a few code examples.  It isn't possible to
say much from just a few formulas.

    Thorsten

____________________________________________________
Thorsten Froehlich, Duisburg, Germany
e-mail: tho### [at] trfde

Visit POV-Ray on the web: http://mac.povray.org

Post a reply to this message

Thorsten Froehlich wrote:

> In article <3DB### [at] freefr> , Le Forgeron <jgr### [at] freefr>  
> wrote:
> 
> 
>>I'm just wondering if any people of the Pov-team reads this newsgroup
>>(and will take the solution from here) or if you need to make a formal
>>posting somewhere else (p.bugreport is moderated, but now that you have
>>discussed here, the moderator should agree to a post... otherwise, I'm
>>afraid of a post in p.general, just to get the attention of the Pov-team
>>members...)
>>
> 
> No patch has been posted.  Just a few code examples.  It isn't possible to
> say much from just a few formulas.
> 


I hope Massimo Valentini will be ok with the following...

Here your requested patch:
      http://jgrimbert.free.fr/pov/correct/blob.pat

I only made it for 3.1, but it should apply fine to 3.5 too!
(once you give the blob.cpp filename to patch when it complains about 
not founding blob.c )

Post a reply to this message