Skip to content

Commit 55f11b0

Browse files
wfpokornywfpokorny
committed
Updating blob accuracy so better synchronized with other shapes.
Fix for GitHub issue POV-Ray#187 Returned minimum intersection depths now set to MIN_ISECT_DEPTH. Internal determine_influences() intersect_<sub-element> functions using >0.0 as was already the subsurface feature. Anything more opens up gaps or jumps where the sub-element density influences are added too late or dropped too earlier. The change does create ray equations slighly more difficult to solve which was likely the reason for the original, largish, v1.0 1e-2 value. The sturm / polysolve solver is sometimes necessary over solve_quartic().
1 parent 27ea23c commit 55f11b0

File tree

2 files changed

+74
-43
lines changed

2 files changed

+74
-43
lines changed

‎source/core/shape/blob.cpp‎

Lines changed: 59 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -119,11 +119,34 @@ namespace pov
119119
* Local preprocessor defines
120120
******************************************************************************/
121121

122-
/* Minimal intersection depth for a valid intersection. */
123-
const DBL DEPTH_TOLERANCE = 1.0e-2;
122+
/// @def DEPTH_TOLERANCE
123+
/// Minimum returned intersection depth for a valid intersection.
124+
///
125+
/// @note
126+
/// Value for @ref DEPTH_TOLERANCE was since v1.0 set to 1.0e-2. Further the
127+
/// value was used for both ling used for the minimum returned and as a mindist
128+
/// in all of the determine_influences()'s intersect_<thing> functions. The
129+
/// latter probably done initially to help the initial solve_quartic() function
130+
/// though using such a large value to determine when to add and remove
131+
/// sub-element influences caused artifacts. In 3.7 when the subsurface (SSLT)
132+
/// feature was added a special conditional existed which set both values to
133+
/// 0.0. Currently a 0.0 value is used for all but the returned intersection
134+
/// where @ref MIN_ISECT_DEPTH is now used. In testing all worked with 0.0, but
135+
/// there is a performance hit to returning all >0.0 intersections. For example,
136+
/// internal refelections at distances which must be ignored in downstream code.
137+
///
138+
const DBL DEPTH_TOLERANCE = MIN_ISECT_DEPTH;
124139

125-
/* Tolerance for inside test. */
126-
const DBL INSIDE_TOLERANCE = 1.0e-6;
140+
/// @def INSIDE_TOLERANCE
141+
/// Density tolerance for inside test.
142+
///
143+
/// @note
144+
/// Value for @ref INSIDE_TOLERANCE was since v1.0 set to 1.0e-6 which is a good
145+
/// value if running single floats. With double floats a value of 1e-15 is better.
146+
/// This is an offset to the inside of the 0.0 density surface after the blob
147+
/// threshold is subtracted. Now using std::numeric_limits to set.
148+
///
149+
const DBL INSIDE_TOLERANCE = 1/std::numeric_limits<DBL>::digits10;
127150

128151
/* Ray enters/exits a component. */
129152
const int ENTERING = 0;
@@ -239,7 +262,6 @@ bool Blob::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadDat
239262
DBL newcoeffs[5], dk[5];
240263
DBL max_bound;
241264
DBL l, w;
242-
DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);
243265

244266
Thread->Stats()[Ray_Blob_Tests]++;
245267

@@ -263,7 +285,7 @@ bool Blob::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadDat
263285

264286
/* Get the intervals along the ray where each component has an effect. */
265287
Blob_Interval_Struct* intervals = Thread->Blob_Intervals;
266-
if ((cnt = determine_influences(P, D, depthTolerance, intervals, Thread)) == 0)
288+
if ((cnt = determine_influences(P, D, intervals, Thread)) == 0)
267289
{
268290
/* Ray doesn't hit any of the component elements. */
269291
return (false);
@@ -528,7 +550,7 @@ bool Blob::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadDat
528550

529551
dist = (dist * max_bound + start_dist) / len;
530552

531-
if ((dist > depthTolerance) && (dist < MAX_DISTANCE))
553+
if ((dist > DEPTH_TOLERANCE) && (dist < MAX_DISTANCE))
532554
{
533555
IPoint = ray.Evaluate(dist);
534556

@@ -600,7 +622,8 @@ bool Blob::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadDat
600622
*
601623
******************************************************************************/
602624

603-
void Blob::insert_hit(const Blob_Element *Element, DBL t0, DBL t1, Blob_Interval_Struct *intervals, unsigned int *cnt)
625+
void Blob::insert_hit(const Blob_Element *Element, DBL t0, DBL t1,
626+
Blob_Interval_Struct *intervals, unsigned int *cnt)
604627
{
605628
unsigned int k;
606629

@@ -678,7 +701,6 @@ void Blob::insert_hit(const Blob_Element *Element, DBL t0, DBL t1, Blob_Interval
678701
*
679702
* Element - Pointer to element structure
680703
* P, D - Ray = P + t * D
681-
* mindist - Min. valid distance
682704
*
683705
* OUTPUT
684706
*
@@ -700,7 +722,8 @@ void Blob::insert_hit(const Blob_Element *Element, DBL t0, DBL t1, Blob_Interval
700722
*
701723
******************************************************************************/
702724

703-
int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
725+
int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P,
726+
const Vector3d& D, DBL *tmin, DBL *tmax)
704727
{
705728
DBL a, b, c, d, t, u, v, w, len;
706729
Vector3d PP, DD;
@@ -786,8 +809,8 @@ int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P, con
786809
*tmin /= len;
787810
*tmax /= len;
788811

789-
if (*tmin < mindist) { *tmin = 0.0; }
790-
if (*tmax < mindist) { *tmax = 0.0; }
812+
if (*tmin < 0.0) { *tmin = 0.0; }
813+
if (*tmax < 0.0) { *tmax = 0.0; }
791814

792815
if (*tmin >= *tmax)
793816
{
@@ -809,7 +832,6 @@ int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P, con
809832
*
810833
* Element - Pointer to element structure
811834
* P, D - Ray = P + t * D
812-
* mindist - Min. valid distance
813835
*
814836
* OUTPUT
815837
*
@@ -831,7 +853,8 @@ int Blob::intersect_cylinder(const Blob_Element *Element, const Vector3d& P, con
831853
*
832854
******************************************************************************/
833855

834-
int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
856+
int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P,
857+
const Vector3d& D, DBL *tmin, DBL *tmax)
835858
{
836859
DBL b, d, t, len;
837860
Vector3d V1, PP, DD;
@@ -855,8 +878,8 @@ int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, co
855878

856879
d = sqrt(d);
857880

858-
*tmax = ( - b + d) / len; if (*tmax < mindist) { *tmax = 0.0; }
859-
*tmin = ( - b - d) / len; if (*tmin < mindist) { *tmin = 0.0; }
881+
*tmax = ( - b + d) / len; if (*tmax < 0.0) { *tmax = 0.0; }
882+
*tmin = ( - b - d) / len; if (*tmin < 0.0) { *tmin = 0.0; }
860883

861884
if (*tmax == *tmin)
862885
{
@@ -885,7 +908,6 @@ int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, co
885908
*
886909
* Element - Pointer to element structure
887910
* P, D - Ray = P + t * D
888-
* mindist - Min. valid distance
889911
*
890912
* OUTPUT
891913
*
@@ -907,7 +929,8 @@ int Blob::intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, co
907929
*
908930
******************************************************************************/
909931

910-
int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
932+
int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P,
933+
const Vector3d& D, DBL *tmin, DBL *tmax)
911934
{
912935
DBL b, d, t, z1, z2, len;
913936
Vector3d PP, DD;
@@ -972,13 +995,13 @@ int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, c
972995
{
973996
/* Ray is crossing the plane from inside to outside. */
974997

975-
*tmin = (t < mindist) ? 0.0 : t;
998+
*tmin = (t < 0.0) ? 0.0 : t;
976999
}
9771000
else
9781001
{
9791002
/* Ray is crossing the plane from outside to inside. */
9801003

981-
*tmax = (t < mindist) ? 0.0 : t;
1004+
*tmax = (t < 0.0) ? 0.0 : t;
9821005
}
9831006

9841007
*tmin /= len;
@@ -1040,13 +1063,13 @@ int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, c
10401063
{
10411064
/* Ray is crossing the plane from inside to outside. */
10421065

1043-
*tmin = (t < mindist) ? 0.0 : t;
1066+
*tmin = (t < 0.0) ? 0.0 : t;
10441067
}
10451068
else
10461069
{
10471070
/* Ray is crossing the plane from outside to inside. */
10481071

1049-
*tmax = (t < mindist) ? 0.0 : t;
1072+
*tmax = (t < 0.0) ? 0.0 : t;
10501073
}
10511074

10521075
*tmin /= len;
@@ -1068,7 +1091,6 @@ int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, c
10681091
*
10691092
* Element - Pointer to element structure
10701093
* P, D - Ray = P + t * D
1071-
* mindist - Min. valid distance
10721094
*
10731095
* OUTPUT
10741096
*
@@ -1090,7 +1112,8 @@ int Blob::intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, c
10901112
*
10911113
******************************************************************************/
10921114

1093-
int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax)
1115+
int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P,
1116+
const Vector3d& D, DBL *tmin, DBL *tmax)
10941117
{
10951118
DBL b, d, t;
10961119
Vector3d V1;
@@ -1108,8 +1131,8 @@ int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P, const
11081131

11091132
d = sqrt(d);
11101133

1111-
*tmax = - b + d; if (*tmax < mindist) { *tmax = 0.0; }
1112-
*tmin = - b - d; if (*tmin < mindist) { *tmin = 0.0; }
1134+
*tmax = - b + d; if (*tmax < 0.0) { *tmax = 0.0; }
1135+
*tmin = - b - d; if (*tmin < 0.0) { *tmin = 0.0; }
11131136

11141137
if (*tmax == *tmin)
11151138
{
@@ -1138,7 +1161,6 @@ int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P, const
11381161
*
11391162
* P, D - Ray = P + t * D
11401163
* Element - Pointer to element structure
1141-
* mindist - Min. valid distance
11421164
*
11431165
* OUTPUT
11441166
*
@@ -1160,7 +1182,8 @@ int Blob::intersect_sphere(const Blob_Element *Element, const Vector3d& P, const
11601182
*
11611183
******************************************************************************/
11621184

1163-
int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Element *Element, DBL mindist, DBL *tmin, DBL *tmax, TraceThreadData *Thread)
1185+
int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Element *Element,
1186+
DBL *tmin, DBL *tmax, TraceThreadData *Thread)
11641187
{
11651188
#ifdef BLOB_EXTRA_STATS
11661189
Thread->Stats()[Blob_Element_Tests]++;
@@ -1173,7 +1196,7 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
11731196
{
11741197
case BLOB_SPHERE:
11751198

1176-
if (!intersect_sphere(Element, P, D, mindist, tmin, tmax))
1199+
if (!intersect_sphere(Element, P, D, tmin, tmax))
11771200
{
11781201
return (false);
11791202
}
@@ -1182,7 +1205,7 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
11821205

11831206
case BLOB_ELLIPSOID:
11841207

1185-
if (!intersect_ellipsoid(Element, P, D, mindist, tmin, tmax))
1208+
if (!intersect_ellipsoid(Element, P, D, tmin, tmax))
11861209
{
11871210
return (false);
11881211
}
@@ -1192,7 +1215,7 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
11921215
case BLOB_BASE_HEMISPHERE:
11931216
case BLOB_APEX_HEMISPHERE:
11941217

1195-
if (!intersect_hemisphere(Element, P, D, mindist, tmin, tmax))
1218+
if (!intersect_hemisphere(Element, P, D, tmin, tmax))
11961219
{
11971220
return (false);
11981221
}
@@ -1201,7 +1224,7 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
12011224

12021225
case BLOB_CYLINDER:
12031226

1204-
if (!intersect_cylinder(Element, P, D, mindist, tmin, tmax))
1227+
if (!intersect_cylinder(Element, P, D, tmin, tmax))
12051228
{
12061229
return (false);
12071230
}
@@ -1228,7 +1251,6 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
12281251
*
12291252
* P, D - Ray = P + t * D
12301253
* Blob - Pointer to blob structure
1231-
* mindist - Min. valid distance
12321254
*
12331255
* OUTPUT
12341256
*
@@ -1253,7 +1275,8 @@ int Blob::intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Ele
12531275
*
12541276
******************************************************************************/
12551277

1256-
int Blob::determine_influences(const Vector3d& P, const Vector3d& D, DBL mindist, Blob_Interval_Struct *intervals, TraceThreadData *Thread) const
1278+
int Blob::determine_influences(const Vector3d& P, const Vector3d& D,
1279+
Blob_Interval_Struct *intervals, TraceThreadData *Thread) const
12571280
{
12581281
unsigned int cnt, size;
12591282
DBL b, t, t0, t1;
@@ -1269,7 +1292,7 @@ int Blob::determine_influences(const Vector3d& P, const Vector3d& D, DBL mindist
12691292

12701293
for (vector<Blob_Element>::iterator i = Data->Entry.begin(); i != Data->Entry.end(); ++i)
12711294
{
1272-
if (intersect_element(P, D, &(*i), mindist, &t0, &t1, Thread))
1295+
if (intersect_element(P, D, &(*i), &t0, &t1, Thread))
12731296
{
12741297
insert_hit(&(*i), t0, t1, intervals, &cnt);
12751298
}
@@ -1293,7 +1316,7 @@ int Blob::determine_influences(const Vector3d& P, const Vector3d& D, DBL mindist
12931316
{
12941317
/* Test element. */
12951318

1296-
if (intersect_element(P, D, reinterpret_cast<Blob_Element *>(Tree->Node), mindist, &t0, &t1, Thread))
1319+
if (intersect_element(P, D, reinterpret_cast<Blob_Element *>(Tree->Node), &t0, &t1, Thread))
12971320
{
12981321
insert_hit(reinterpret_cast<Blob_Element *>(Tree->Node), t0, t1, intervals, &cnt);
12991322
}

‎source/core/shape/blob.h‎

Lines changed: 15 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -170,16 +170,24 @@ class Blob : public ObjectBase
170170
static void Transform_Blob_Element(Blob_Element *Element, const TRANSFORM *Trans);
171171
private:
172172
static void element_normal(Vector3d& Result, const Vector3d& P, const Blob_Element *Element);
173-
static int intersect_element(const Vector3d& P, const Vector3d& D, const Blob_Element *Element, DBL mindist, DBL *t0, DBL *t1, TraceThreadData *Thread);
174-
static void insert_hit(const Blob_Element *Element, DBL t0, DBL t1, Blob_Interval_Struct *intervals, unsigned int *cnt);
175-
int determine_influences(const Vector3d& P, const Vector3d& D, DBL mindist, Blob_Interval_Struct *intervals, TraceThreadData *Thread) const;
173+
static int intersect_element(const Vector3d& P, const Vector3d& D,
174+
const Blob_Element *Element, DBL *t0, DBL *t1,
175+
TraceThreadData *Thread);
176+
static void insert_hit(const Blob_Element *Element, DBL t0, DBL t1,
177+
Blob_Interval_Struct *intervals, unsigned int *cnt);
178+
int determine_influences(const Vector3d& P, const Vector3d& D,
179+
Blob_Interval_Struct *intervals, TraceThreadData *Thread) const;
176180
DBL calculate_field_value(const Vector3d& P, TraceThreadData *Thread) const;
177181
static DBL calculate_element_field(const Blob_Element *Element, const Vector3d& P);
178182

179-
static int intersect_cylinder(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax);
180-
static int intersect_hemisphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax);
181-
static int intersect_sphere(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax);
182-
static int intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P, const Vector3d& D, DBL mindist, DBL *tmin, DBL *tmax);
183+
static int intersect_cylinder(const Blob_Element *Element, const Vector3d& P,
184+
const Vector3d& D, DBL *tmin, DBL *tmax);
185+
static int intersect_hemisphere(const Blob_Element *Element, const Vector3d& P,
186+
const Vector3d& D, DBL *tmin, DBL *tmax);
187+
static int intersect_sphere(const Blob_Element *Element, const Vector3d& P,
188+
const Vector3d& D, DBL *tmin, DBL *tmax);
189+
static int intersect_ellipsoid(const Blob_Element *Element, const Vector3d& P,
190+
const Vector3d& D, DBL *tmin, DBL *tmax);
183191

184192
static void get_element_bounding_sphere(const Blob_Element *Element, Vector3d& Center, DBL *Radius2);
185193
void build_bounding_hierarchy();

0 commit comments

Comments
 (0)