I faced this bugreport and hope to have found where's the problem.

It seems that when computing the normal of the superellipsoid some
kind of over/underflow occurs with n or e very small, reordering the
computation produces a better image. Here's the diff to be applied
to super.cpp and a scene that shows the difference.

Massimo

#macro S(E,N,P)
  superellipsoid { <E, N> translate P scale 2*y+1 pigment {rgb
<1,1,0>} }
#end
light_source { -35*z, rgb 1 }
camera {location -35*z look_at 0 }

#declare E = 5;
#declare P = -20*x;
#while (E > 0.00001)
  S(E,.1,P+3.5*y)
  S(E, E, P)
  S(.1,E,P-3.5*y)
  #declare P = P + 5*x;
  #declare E = E / 5;
#end


479,480d478
<   DBL k, x, y, z;
<   VECTOR P, N, E;
481a480,481
>   VECTOR const& E = Superellipsoid->Power;
>   VECTOR P;
487,498c487,494
<   Assign_Vector(E, Superellipsoid->Power);
<
<   x = fabs(P[X]);
<   y = fabs(P[Y]);
<   z = fabs(P[Z]);
<
<   k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
<
<   N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
<   N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
<   N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
<
---
>   DBL x = P[X] != 0 ? power(fabs(P[X]), E[X]) : 0;
>   DBL y = P[Y] != 0 ? power(fabs(P[Y]), E[X]) : 0;
>   DBL z = P[Z] != 0 ? power(fabs(P[Z]), E[Z]) : 0;
>
>   P[X] = x ? (1 - z) * x / P[X] : 0;
>   P[Y] = y ? (1 - z) * y / P[Y] : 0;
>   P[Z] = z ? (x+y) ? (x + y) * z / P[Z] : 1 : 0;
>
501c497
<   MTransNormal(Result, N, Superellipsoid->Trans);
---
>   MTransNormal(Result, P, Superellipsoid->Trans);

Post a reply to this message

"Massimo Valentini" ha scritto 
: 
: #macro S(E,N,P)
:   superellipsoid { <E, N> translate P scale 2*y+1 pigment {rgb
: <1,1,0>} }
: #end
: light_source { -35*z, rgb 1 }
: camera {location -35*z look_at 0 }
: 
: #declare E = 5;
: #declare P = -20*x;
: #while (E > 0.00001)
:   S(E,.1,P+3.5*y)
:   S(E, E, P)
:   S(.1,E,P-3.5*y)
:   #declare P = P + 5*x;
:   #declare E = E / 5;
: #end
: 

With the following patch to the original file super.cpp it's
possible to obtain a better image for the previous scene.

I modified the even the function evaluate_superellipsoid too and
introduced a helper function.

HTH Massimo


479,480d478
<   DBL k, x, y, z;
<   VECTOR P, N, E;
481a480,481
>   VECTOR const& E = Superellipsoid->Power;
>   VECTOR P, N;
487,491c487,492
<   Assign_Vector(E, Superellipsoid->Power);
< 
<   x = fabs(P[X]);
<   y = fabs(P[Y]);
<   z = fabs(P[Z]);
---
>   DBL r, z2n = 0;
>   if (P[Z] != 0)
>   {
>     z2n = power(fabs(P[Z]), E[Z]);
>     P[Z] = z2n / P[Z];
>   }
493c494,496
<   k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
---
>   if (fabs(P[X]) > fabs(P[Y]))
>   {
>     r = power(fabs(P[Y] / P[X]), E[X]);
495,497c498,503
<   N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
<   N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
<   N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
---
>     P[X] = (1-z2n)  /  P[X];
>     P[Y] = P[Y] ? (1-z2n) * r / P[Y] : 0;
>   }
>   else if (P[Y] != 0)
>   {
>     r = power(fabs(P[X] / P[Y]), E[X]);
498a505,508
>     P[X] = P[X] ? (1-z2n) * r / P[X] : 0;
>     P[Y] = (1-z2n) / P[Y];
>   }
>   if (P[Z]) P[Z] *= (1 + r);
501c511
<   MTransNormal(Result, N, Superellipsoid->Trans);
---
>   MTransNormal(Result, P, Superellipsoid->Trans);
1084a1095,1111
> static DBL my_g(DBL x, DBL y, DBL e)
> {
>   DBL g = 0;
>   if (x > y)
>   {
>     g = 1 + power(y/x, e);
>     if (g != 1) g = power(g, 1/e);
>     g *= x;
>   }
>   else if (y != 0)
>   {
>     g = 1 + power(x/y, e);
>     if (g != 1) g = power(g, 1/e);
>     g *= y;
>   }
>   return g;
> }
1088,1094c1115,1117
<   VECTOR V1;
< 
<   V1[X] = power(fabs(P[X]), Superellipsoid->Power[X]);
<   V1[Y] = power(fabs(P[Y]), Superellipsoid->Power[X]);
<   V1[Z] = power(fabs(P[Z]), Superellipsoid->Power[Z]);
< 
<   return(power(V1[X] + V1[Y], Superellipsoid->Power[Y]) + V1[Z] - 1.0);
---
>   DBL f = my_g(fabs(P[X]), fabs(P[Y]), Superellipsoid->Power[X]);
>   f = my_g(f, fabs(P[Z]), Superellipsoid->Power[Z]);
>   return f - 1;

Post a reply to this message

On Thu, 17 Oct 2002 12:02:11 +0200, "Massimo Valentini"
<six### [at] inwindit> wrote:
>479,480d478
><   DBL k, x, y, z;
><   VECTOR P, N, E;

lots of stuff snipped.

If you want others to be able to use your diff you better make it in
context format, 3 lines before/after changes used to be the rule.
Otherwise your diff can only be patched to the exact same source level
you diff'd with.

Back in the old band-width constrained days of POV-Ray 1.0 development
we often went thru 10-20 patch levels between full source syncs and
that absolutely requires context diffs. As a side note the FracInt
boys often went to 50+ patch levels between full source syncs, but
that's a bit extreme.

/Erkki

Post a reply to this message

:
: If you want others to be able to use your diff you better make it in
: context format, 3 lines before/after changes used to be the rule.
: Otherwise your diff can only be patched to the exact same source level
: you diff'd with.
:
Here is a diff in the unified format (3 lines before/after changes).
Massimo

--- povray-3.50a/src/super.cpp Fri Oct 18 12:45:32 2002
+++ povray-3.50b/src/super.cpp Tue Oct 22 11:18:09 2002
@@ -476,29 +476,39 @@

 static void Superellipsoid_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
 {
-  DBL k, x, y, z;
-  VECTOR P, N, E;
   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
+  VECTOR const& E = Superellipsoid->Power;
+  VECTOR P, N;

   /* Transform the point into the superellipsoid space. */

   MInvTransPoint(P, Inter->IPoint, Superellipsoid->Trans);

-  Assign_Vector(E, Superellipsoid->Power);
-
-  x = fabs(P[X]);
-  y = fabs(P[Y]);
-  z = fabs(P[Z]);
+  DBL r, z2n = 0;
+  if (P[Z] != 0)
+  {
+    z2n = power(fabs(P[Z]), E[Z]);
+    P[Z] = z2n / P[Z];
+  }

-  k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
+  if (fabs(P[X]) > fabs(P[Y]))
+  {
+    r = power(fabs(P[Y] / P[X]), E[X]);

-  N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
-  N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
-  N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
+    P[X] = (1-z2n)  /  P[X];
+    P[Y] = P[Y] ? (1-z2n) * r / P[Y] : 0;
+  }
+  else if (P[Y] != 0)
+  {
+    r = power(fabs(P[X] / P[Y]), E[X]);

+    P[X] = P[X] ? (1-z2n) * r / P[X] : 0;
+    P[Y] = (1-z2n) / P[Y];
+  }
+  if (P[Z]) P[Z] *= (1 + r);
   /* Transform the normalt out of the superellipsoid space. */

-  MTransNormal(Result, N, Superellipsoid->Trans);
+  MTransNormal(Result, P, Superellipsoid->Trans);

   VNormalize(Result, Result);
 }
@@ -1082,16 +1092,29 @@
 *   Oct 1994 : Creation.
 *
 ******************************************************************************/
+static DBL my_g(DBL x, DBL y, DBL e)
+{
+  DBL g = 0;
+  if (x > y)
+  {
+    g = 1 + power(y/x, e);
+    if (g != 1) g = power(g, 1/e);
+    g *= x;
+  }
+  else if (y != 0)
+  {
+    g = 1 + power(x/y, e);
+    if (g != 1) g = power(g, 1/e);
+    g *= y;
+  }
+  return g;
+}

 static DBL evaluate_superellipsoid(VECTOR P, SUPERELLIPSOID *Superellipsoid)
 {
-  VECTOR V1;
-
-  V1[X] = power(fabs(P[X]), Superellipsoid->Power[X]);
-  V1[Y] = power(fabs(P[Y]), Superellipsoid->Power[X]);
-  V1[Z] = power(fabs(P[Z]), Superellipsoid->Power[Z]);
-
-  return(power(V1[X] + V1[Y], Superellipsoid->Power[Y]) + V1[Z] - 1.0);
+  DBL f = my_g(fabs(P[X]), fabs(P[Y]), Superellipsoid->Power[X]);
+  f = my_g(f, fabs(P[Z]), Superellipsoid->Power[Z]);
+  return f - 1;
 }

Post a reply to this message