Fix bulk of GitHub issue #147 sphere sweep artifacts for good.

Author: c-lipka

changes.txt

@@ -13,6 +13,11 @@ Notes:
   For example, to read FS#270, visit:
   http://bugs.povray.org/270
 
+  (Note: As of 3.7.1-beta.8 our flyspray bug tracker is officially mothballed,
+  and we no longer refer to it in the changelog. Any flyspray reports that
+  were still open at that point have been duplicated as GitHub issues titled
+  "Port of FSxxx - ...")
+
 - To view a referenced newsgroup posting (<[email protected]>), prefix
   the message ID with the URL http://news.povray.org/.
   The '<' and '>' are optional (if using a shell you may want to omit them).
@@ -27,7 +32,15 @@ should be located in the same directory as this file is.
 Changes between 3.7.1-beta.7 and 3.7.1-TBD
 ==========================================
 
-Miscellaneous:
+Fixed or Mitigated Bugs
+-----------------------
+
+Reported via GitHub:
+
+  - #147 ("Port of FS81 - sphere_sweep generating artifacts")
+
+Other Noteworthy
+----------------
 
   - Minor updates to the documentation.
 

revision.txt

@@ -19,13 +19,22 @@ Notes:
   For example, to read FS#270, visit:
   http://bugs.povray.org/270
 
+  (Note: As of 3.7.1-beta.8 our flyspray bug tracker is officially mothballed,
+  and we no longer refer to it in the changelog. Any flyspray reports that
+  were still open at that point have been duplicated as GitHub issues titled
+  "Port of FSxxx - ...")
+
 - To view a referenced newsgroup posting (<[email protected]>), prefix
   the message ID with the URL http://news.povray.org/.
   The '<' and '>' are optional (if using a shell you may want to omit them).
   For example, to read <[email protected]>, visit:
   http://news.povray.org/[email protected]
 ------------------------------------------------------------------------------
 
+------------------------------------------------------------------------------
+POV-Ray v3.7.1-beta.TBD
+------------------------------------------------------------------------------
+
 ------------------------------------------------------------------------------
 POV-Ray v3.7.1-beta.7
 ------------------------------------------------------------------------------

source/core/shape/spheresweep.cpp

@@ -101,6 +101,7 @@ namespace pov
 
 const DBL DEPTH_TOLERANCE       = 1.0e-6;
 const DBL ZERO_TOLERANCE        = 1.0e-4;
+const DBL RADIUS_TOLERANCE      = 1.0e-6;
 
 // TODO - Technically this should be `(Num_Spheres + Num_Segments * order) * 2` [CLi]
 #define SPHSWEEP_MAX_ISECT  ((Num_Spheres + Num_Segments) * 2 + 64)
@@ -408,7 +409,6 @@ int SphereSweep::Intersect_Segment(const BasicRay &ray, const SPHSWEEP_SEG *Segm
     DBL             fp0, fp1;
     DBL             t;
     SPHSWEEP_SPH    Temp_Sphere;
-    Vector3d        Center;
 
     Isect_Count = 0;
 
@@ -617,51 +617,31 @@ int SphereSweep::Intersect_Segment(const BasicRay &ray, const SPHSWEEP_SEG *Segm
             m++;
     }
 
-    switch(Segment->Num_Coefs)
+    for (m = 0; m < Num_Poly_Roots; m++)
     {
-        case 2:
-            for(m = 0; m < Num_Poly_Roots; m++)
-            {
-                fp0 = 2.0 * d * Root[m]
-                    + e;
-                fp1 = b;
-
-                if(fabs(fp1) > ZERO_TOLERANCE)
-                {
-                    t = -fp0 / fp1;
-
-                    if((t > -MAX_DISTANCE) && (t < MAX_DISTANCE))
-                    {
-                        // Add intersection
-                        Isect[Isect_Count].t = t;
+        // Calculate center and radius of sphere at root
 
-                        // Calculate point
-                        Isect[Isect_Count].Point = ray.Evaluate(t);
+        DBL rootPwr = 1;
+        Temp_Sphere.Center = Segment->Center_Coef[0];
+        Temp_Sphere.Radius = Segment->Radius_Coef[0];
+        for (n = 1; n < Segment->Num_Coefs; ++n)
+        {
+            rootPwr *= Root[m];
+            Temp_Sphere.Center += Segment->Center_Coef[n] * rootPwr;
+            Temp_Sphere.Radius += Segment->Radius_Coef[n] * rootPwr;
+        }
 
-                        // Calculate normal
-                        Center = Segment->Center_Coef[0] + Root[m] * Segment->Center_Coef[1];
-                        Isect[Isect_Count].Normal = (Isect[Isect_Count].Point - Center).normalized();
+        switch (Segment->Num_Coefs)
+        {
+            case 2:
 
-                        Isect_Count++;
-                    }
-                }
-                else
-                {
-                    // Calculate center of single sphere
-                    Temp_Sphere.Center = Segment->Center_Coef[1] * Root[m] + Segment->Center_Coef[0];
+                fp0 = 2.0 * d * Root[m]
+                    + e;
+                fp1 = b;
+                break;
 
-                    // Calculate radius of single sphere
-                    Temp_Sphere.Radius = Segment->Radius_Coef[1] * Root[m] + Segment->Radius_Coef[0];
+            case 4:
 
-                    // Calculate intersections
-                    if(Intersect_Sphere(ray, &Temp_Sphere, Isect + Isect_Count))
-                        Isect_Count += 2;
-                }
-            }
-            break;
-        case 4:
-            for (m = 0; m < Num_Poly_Roots; m++)
-            {
                 fp0 = 6.0 * f * Root[m] * Root[m] * Root[m] * Root[m] * Root[m]
                     + 5.0 * g * Root[m] * Root[m] * Root[m] * Root[m]
                     + 4.0 * h * Root[m] * Root[m] * Root[m]
@@ -671,49 +651,49 @@ int SphereSweep::Intersect_Segment(const BasicRay &ray, const SPHSWEEP_SEG *Segm
                 fp1 = 3.0 * b * Root[m] * Root[m]
                     + 2.0 * c * Root[m]
                     + d;
+                break;
 
-                if(fabs(fp1) > ZERO_TOLERANCE)
-                {
-                    t = -fp0 / fp1;
+            default:
+                POV_SHAPE_ASSERT(false);
+                break;
+        }
 
-                    if((t > -MAX_DISTANCE) && (t < MAX_DISTANCE))
-                    {
-                        // Add intersection
-                        Isect[Isect_Count].t = t;
+        bool doubleIsect = (fabs(fp1) < ZERO_TOLERANCE);
+        if (!doubleIsect)
+        {
+            // Looks like a single-intersection root at first glance
 
-                        // Calculate point
-                        Isect[Isect_Count].Point = ray.Evaluate(t);
+            t = -fp0 / fp1;
 
-                        // Calculate normal
-                        Center = Segment->Center_Coef[3] * (Root[m] * Root[m] * Root[m])
-                               + Segment->Center_Coef[2] * (Root[m] * Root[m])
-                               + Segment->Center_Coef[1] *  Root[m]
-                               + Segment->Center_Coef[0];
-                        Isect[Isect_Count].Normal = (Isect[Isect_Count].Point - Center).normalized();
+            if ((t > -MAX_DISTANCE) && (t < MAX_DISTANCE))
+            {
+                // Calculate intersection point
+                Isect[Isect_Count].Point = ray.Evaluate(t);
 
-                        Isect_Count++;
-                    }
-                }
-                else
+                // Calculate normal
+                Isect[Isect_Count].Normal = (Isect[Isect_Count].Point - Temp_Sphere.Center).normalized();
+
+                // Sanity-check intersection point for proper radius
+                DBL rIsect = (Isect[Isect_Count].Point - Temp_Sphere.Center).length();
+                DBL rRoot  = Temp_Sphere.Radius;
+
+                // Check for numeric instability; this may (only?) happen if the root corresponds to two intersections
+                doubleIsect = fabs(rIsect - rRoot) > RADIUS_TOLERANCE;
+                if (!doubleIsect)
                 {
-                    // Calculate center of single sphere
-                    Temp_Sphere.Center = Segment->Center_Coef[3] * (Root[m] * Root[m] * Root[m])
-                                       + Segment->Center_Coef[2] * (Root[m] * Root[m])
-                                       + Segment->Center_Coef[1] *  Root[m]
-                                       + Segment->Center_Coef[0];
-
-                    // Calculate radius of single sphere
-                    Temp_Sphere.Radius = Segment->Radius_Coef[3] * Root[m] * Root[m] * Root[m]
-                                       + Segment->Radius_Coef[2] * Root[m] * Root[m]
-                                       + Segment->Radius_Coef[1] * Root[m]
-                                       + Segment->Radius_Coef[0];
-
-                    // Calculate intersections
-                    if (Intersect_Sphere(ray, &Temp_Sphere, Isect + Isect_Count))
-                        Isect_Count += 2;
+                    // Add intersection
+                    Isect[Isect_Count].t = t;
+                    Isect_Count++;
                 }
             }
-            break;
+        }
+
+        if (doubleIsect)
+        {
+            // Calculate intersections
+            if(Intersect_Sphere(ray, &Temp_Sphere, Isect + Isect_Count))
+                Isect_Count += 2;
+        }
     }
 
     return Isect_Count;
@@ -1353,10 +1333,8 @@ void SphereSweep::Compute_BBox()
     {
         for (i = 0; i < Num_Modeling_Spheres; i++)
         {
-            {
-                mins = min(mins, Modeling_Sphere[i].Center - Vector3d(fabs(Modeling_Sphere[i].Radius)));
-                maxs = max(maxs, Modeling_Sphere[i].Center + Vector3d(fabs(Modeling_Sphere[i].Radius)));
-            }
+            mins = min(mins, Modeling_Sphere[i].Center - Vector3d(fabs(Modeling_Sphere[i].Radius)));
+            maxs = max(maxs, Modeling_Sphere[i].Center + Vector3d(fabs(Modeling_Sphere[i].Radius)));
         }
     }