Implemented potential pattern and isosurface polarity keyword. (#88)

Author: c-lipka

changes.txt

@@ -49,6 +49,12 @@ Prior to the release of 3.7.1, the following items still need urgent attention:
 New Features
 ------------
 
+- A new pattern, `potential`, has been added to define a pattern based on the
+  potential field of a blob or isosurface object.
+
+- A new parameter, `polarity`, has been added to isosurfaces to choose whether
+  function values above the threshold should indicate the inside or outside.
+
 - Basic support for importing Wavefront OBJ files has been added as an extension
   to the `mesh` primitive.
 

source/base/version.h

@@ -45,7 +45,7 @@
 #define OFFICIAL_VERSION_STRING "3.7.1"
 #define OFFICIAL_VERSION_NUMBER 371
 
-#define POV_RAY_PRERELEASE "x.obj.8787755"
+#define POV_RAY_PRERELEASE "x.potential.8787796"
 
 #if (POV_RAY_IS_AUTOBUILD == 1) && ((POV_RAY_IS_OFFICIAL == 1) || (POV_RAY_IS_SEMI_OFFICIAL == 1))
 #ifdef POV_RAY_PRERELEASE

source/core/material/pattern.cpp

@@ -301,6 +301,9 @@ DBL ContinuousPattern::Evaluate(const Vector3d& EPoint, const Intersection *pIse
 {
     DBL value = EvaluateRaw(EPoint, pIsection, pRay, pThread);
 
+    if (waveType == kWaveType_Raw)
+        return value;
+
     if(waveFrequency != 0.0)
         value = fmod(value * waveFrequency + wavePhase, 1.00001); // TODO FIXME - magic number! Should be 1.0+SOME_EPSILON (or maybe actually 1.0?)
 
@@ -467,6 +470,29 @@ PigmentPattern::~PigmentPattern()
 }
 
 
+PotentialPattern::PotentialPattern() :
+    pObject(NULL),
+    subtractThreshold(false)
+{
+    waveType = kWaveType_Raw;
+}
+
+PotentialPattern::PotentialPattern(const PotentialPattern& obj) :
+    ContinuousPattern(obj),
+    pObject(NULL),
+    subtractThreshold(obj.subtractThreshold)
+{
+    if (obj.pObject)
+        pObject = Copy_Object(obj.pObject);
+}
+
+PotentialPattern::~PotentialPattern()
+{
+    if (pObject)
+        Destroy_Object(pObject);
+}
+
+
 ColourBlendMapConstPtr RadialPattern::GetDefaultBlendMap() const { return gpDefaultBlendMap_Radial; }
 
 ColourBlendMapConstPtr SquarePattern::GetDefaultBlendMap() const { return gpDefaultBlendMap_Square; }
@@ -7802,6 +7828,17 @@ DBL PlanarPattern::EvaluateRaw(const Vector3d& EPoint, const Intersection *pIsec
 }
 
 
+/*****************************************************************************/
+
+DBL PotentialPattern::EvaluateRaw(const Vector3d& EPoint, const Intersection *pIsection, const Ray *pRay, TraceThreadData *pThread) const
+{
+    if (pObject != NULL)
+        return pObject->GetPotential (EPoint, subtractThreshold, pThread);
+
+    return 0.0;
+}
+
+
 /*****************************************************************************
 *
 * FUNCTION

source/core/material/pattern.h

@@ -127,6 +127,7 @@ enum PATTERN_IDS
 
 enum WaveType
 {
+    kWaveType_Raw,      ///< Use raw pattern value, allowing it to exceed the [0..1] range.
     kWaveType_Ramp,     ///< Ramps up from 0 to 1, then drops sharply back to 0, and repeats.
     kWaveType_Sine,     ///< Oscillates between 0 and 1 using a sine-based formula.
     kWaveType_Triangle, ///< Ramps up from 0 to 1, then ramps down from 1 to 0, and repeats.
@@ -693,6 +694,22 @@ struct PlanarPattern : public ContinuousPattern
     virtual DBL EvaluateRaw(const Vector3d& EPoint, const Intersection *pIsection, const Ray *pRay, TraceThreadData *pThread) const;
 };
 
+/// Implements the `potential` pattern.
+///
+/// @todo   The additional member variables should possibly be encapsulated.
+///
+struct PotentialPattern : public ContinuousPattern
+{
+    ObjectPtr   pObject;
+    bool        subtractThreshold;
+
+    PotentialPattern();
+    PotentialPattern(const PotentialPattern& obj);
+    virtual ~PotentialPattern();
+    virtual PatternPtr Clone() const { return BasicPattern::Clone(*this); }
+    virtual DBL EvaluateRaw(const Vector3d& EPoint, const Intersection *pIsection, const Ray *pRay, TraceThreadData *pThread) const;
+};
+
 /// Implements the `quilted` pattern.
 struct QuiltedPattern : public ContinuousPattern
 {

source/core/scene/object.cpp

@@ -826,10 +826,16 @@ void Destroy_Object(ObjectPtr Object)
 
 ObjectBase::~ObjectBase()
 {
-	Destroy_Transform(Trans);
+    Destroy_Transform(Trans);
 }
 
 
+double ObjectBase::GetPotential (const Vector3d& p, bool subtractThreshold, TraceThreadData *threaddata) const
+{
+    POV_ASSERT (false);
+    return 0.0;
+}
+
 /*****************************************************************************
 *
 * FUNCTION

source/core/scene/object.h

@@ -134,23 +134,24 @@ namespace pov
 
 /* Object types. */
 
-#define BASIC_OBJECT                0
-#define PATCH_OBJECT                1 /* Has no inside, no inverse */
-#define TEXTURED_OBJECT             2 /* Has texture, possibly in children */
-#define IS_COMPOUND_OBJECT          4 /* Has children field */
-#define STURM_OK_OBJECT             8 /* STURM legal */
-//#define WATER_LEVEL_OK_OBJECT      16 /* WATER_LEVEL legal */
-#define LIGHT_SOURCE_OBJECT        32 /* link me in frame.light_sources */
-#define BOUNDING_OBJECT            64 /* This is a holder for bounded object */
-//#define SMOOTH_OK_OBJECT          128 /* SMOOTH legal */
-#define IS_CHILD_OBJECT           256 /* Object is inside a COMPOUND */
+#define BASIC_OBJECT                0x0000u
+#define PATCH_OBJECT                0x0001u // Has no inside, no inverse
+#define TEXTURED_OBJECT             0x0002u // Has texture, possibly in children
+#define IS_COMPOUND_OBJECT          0x0004u // Has children field
+#define STURM_OK_OBJECT             0x0008u /* STURM legal */
+//#define WATER_LEVEL_OK_OBJECT         0x0010u /* WATER_LEVEL legal */
+#define LIGHT_SOURCE_OBJECT         0x0020u /* link me in frame.light_sources */
+#define BOUNDING_OBJECT             0x0040u /* This is a holder for bounded object */
+//#define SMOOTH_OK_OBJECT              0x0080u /* SMOOTH legal */
+#define IS_CHILD_OBJECT             0x0100u /* Object is inside a COMPOUND */
 /* NK 1998 - DOUBLE_ILLUMINATE is not used anymore - use DOUBLE_ILLUMINATE_FLAG */
-#define HIERARCHY_OK_OBJECT       512 /* NO_HIERARCHY legal */
-#define LT_SRC_UNION_OBJECT      1024 /* Union of light_source objects only */
-#define LIGHT_GROUP_OBJECT       2048 /* light_group union object [trf] */
-#define LIGHT_GROUP_LIGHT_OBJECT 4096 /* light in light_group object [trf] */
-#define CSG_DIFFERENCE_OBJECT    8192 /* csg difference object */
-#define IS_CSG_OBJECT           16384 /* object is a csg and not some other compound object */
+#define HIERARCHY_OK_OBJECT         0x0200u /* NO_HIERARCHY legal */
+#define LT_SRC_UNION_OBJECT         0x0400u /* Union of light_source objects only */
+#define LIGHT_GROUP_OBJECT          0x0800u /* light_group union object [trf] */
+#define LIGHT_GROUP_LIGHT_OBJECT    0x1000u /* light in light_group object [trf] */
+#define CSG_DIFFERENCE_OBJECT       0x2000u /* csg difference object */
+#define IS_CSG_OBJECT               0x4000u /* object is a csg and not some other compound object */
+#define POTENTIAL_OBJECT            0x8000u ///< Object has an intrinsic potential field associated.
 #define CHILDREN_FLAGS (PATCH_OBJECT+TEXTURED_OBJECT)  /* Reverse inherited flags */
 
 
@@ -222,6 +223,7 @@ class ObjectBase
         virtual bool Precompute() { return true; };
 
         virtual bool All_Intersections(const Ray&, IStack&, TraceThreadData *) = 0; // could be "const", if it wasn't for isosurface max_gradient estimation stuff
+        virtual double GetPotential (const Vector3d&, bool subtractThreshold, TraceThreadData *) const;
         virtual bool Inside(const Vector3d&, TraceThreadData *) const = 0;
         virtual void Normal(Vector3d&, Intersection *, TraceThreadData *) const = 0;
         virtual void UVCoord(Vector2d&, const Intersection *, TraceThreadData *) const;

source/core/shape/blob.cpp

@@ -1610,6 +1610,25 @@ bool Blob::Inside(const Vector3d& Test_Point, TraceThreadData *Thread) const
 }
 
 
+double Blob::GetPotential (const Vector3d& globalPoint, bool subtractThreshold, TraceThreadData *pThread) const
+{
+    Vector3d localPoint;
+
+    if (Trans != NULL)
+        MInvTransPoint (localPoint, globalPoint, Trans);
+    else
+        localPoint = globalPoint;
+
+    double potential = calculate_field_value (localPoint, pThread);
+    if (subtractThreshold)
+        potential -= Data->Threshold;
+
+    if (Test_Flag (this, INVERTED_FLAG))
+        return -potential;
+    else
+        return  potential;
+}
+
 
 /*****************************************************************************
 *

source/core/shape/blob.h

@@ -49,7 +49,7 @@ namespace pov
 * Global preprocessor defines
 ******************************************************************************/
 
-#define BLOB_OBJECT (STURM_OK_OBJECT+HIERARCHY_OK_OBJECT)
+#define BLOB_OBJECT (STURM_OK_OBJECT+HIERARCHY_OK_OBJECT+POTENTIAL_OBJECT)
 
 // TODO - the following values probably don't have to be all discrete bits, except for BLOB_ENTER_EXIT_FLAG
 #define BLOB_ENTER_EXIT_FLAG      1 // internal use only
@@ -136,6 +136,7 @@ class Blob : public ObjectBase
 
         virtual bool All_Intersections(const Ray&, IStack&, TraceThreadData *);
         virtual bool Inside(const Vector3d&, TraceThreadData *) const;
+        virtual double GetPotential (const Vector3d&, bool subtractThreshold, TraceThreadData *) const;
         virtual void Normal(Vector3d&, Intersection *, TraceThreadData *) const;
         virtual void Translate(const Vector3d&, const TRANSFORM *);
         virtual void Rotate(const Vector3d&, const TRANSFORM *);

source/core/shape/isosurface.cpp

@@ -128,7 +128,7 @@ bool IsoSurface::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThr
         if(closed != false)
         {
             VTmp = Plocal + Depth1 * Dlocal;
-            tmp = fn.Evaluate(VTmp) - threshold;
+            tmp = EvaluatePolarized (fn, VTmp);
             if(Depth1 > accuracy)
             {
                 if(tmp < 0.0)                   /* The ray hits the bounding shape */
@@ -149,12 +149,12 @@ bool IsoSurface::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThr
                 {
                     Depth1 = accuracy * 5.0;
                     VTmp = Plocal + Depth1 * Dlocal;
-                    if(fn.Evaluate(VTmp) < threshold)
+                    if (IsInside (fn, VTmp))
                         Thread->isosurfaceData->Inv3 = -1;
                     /* Change the sign of the function (IPoint is in the bounding shpae.)*/
                 }
                 VTmp = Plocal + Depth2 * Dlocal;
-                if(fn.Evaluate(VTmp) < threshold)
+                if (IsInside (fn, VTmp))
                 {
                     IPoint = ray.Evaluate(Depth2);
                     if(Clip.empty() || Point_In_Clip(IPoint, Clip, Thread))
@@ -180,11 +180,11 @@ bool IsoSurface::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThr
         {
             /* IPoint is on the isosurface */
             VTmp = Plocal + tmin * Dlocal;
-            if(fabs(fn.Evaluate(VTmp) - threshold) < (maxg * accuracy * 4.0))
+            if (EvaluateAbs (fn, VTmp) < (maxg * accuracy * 4.0))
             {
                 tmin = accuracy * 5.0;
                 VTmp = Plocal + tmin * Dlocal;
-                if(fn.Evaluate(VTmp) < threshold)
+                if (IsInside (fn, VTmp))
                     Thread->isosurfaceData->Inv3 = -1;
                 /* change the sign and go into the isosurface */
             }
@@ -213,6 +213,8 @@ bool IsoSurface::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThr
 
         if(IFound)
             Thread->Stats()[Ray_IsoSurface_Tests_Succeeded]++;
+
+        Thread->isosurfaceData->pFn = NULL;
     }
 
     if(eval == true)
@@ -264,14 +266,34 @@ bool IsoSurface::Inside(const Vector3d& IPoint, TraceThreadData *Thread) const
     if(!container->Inside(New_Point))
         return (Test_Flag(this, INVERTED_FLAG));
 
-    if(GenericScalarFunctionInstance(Function, Thread).Evaluate(New_Point) > threshold)
+    GenericScalarFunctionInstance fn(Function, Thread);
+    if (!IsInside (fn, New_Point))
         return (Test_Flag(this, INVERTED_FLAG));
 
     /* Inside the box. */
     return (!Test_Flag(this, INVERTED_FLAG));
 }
 
 
+double IsoSurface::GetPotential (const Vector3d& globalPoint, bool subtractThreshold, TraceThreadData *pThread) const
+{
+    Vector3d localPoint;
+
+    if (Trans != NULL)
+        MInvTransPoint (localPoint, globalPoint, Trans);
+    else
+        localPoint = globalPoint;
+
+    double potential = GenericScalarFunctionInstance(Function, pThread).Evaluate (localPoint);
+    if (subtractThreshold)
+        potential -= threshold;
+
+    if (Test_Flag (this, INVERTED_FLAG))
+        return -potential;
+    else
+        return  potential;
+}
+
 
 /*****************************************************************************
 *
@@ -345,6 +367,9 @@ void IsoSurface::Normal(Vector3d& Result, Intersection *Inter, TraceThreadData *
 
             Result.normalize();
         }
+
+        if (positivePolarity)
+            Result.invert();
     }
 }
 
@@ -513,7 +538,9 @@ void IsoSurface::Transform(const TRANSFORM* tr)
 *
 ******************************************************************************/
 
-IsoSurface::IsoSurface() : ObjectBase(ISOSURFACE_OBJECT)
+IsoSurface::IsoSurface() :
+    ObjectBase(ISOSURFACE_OBJECT),
+    positivePolarity(false)
 {
     container = shared_ptr<ContainedByShape>(new ContainedByBox());
 
@@ -588,6 +615,8 @@ ObjectPtr IsoSurface::Copy()
     // mark it as copy for use by max_gradient warning code
     New->isCopy = true;
 
+    New->positivePolarity = positivePolarity;
+
     New->container = shared_ptr<ContainedByShape>(container->Copy());
 
     return (New);
@@ -964,7 +993,31 @@ DBL IsoSurface::Float_Function(ISO_ThreadData& itd, DBL t) const
 
     VTmp = itd.Pglobal + t * itd.Dglobal;
 
-    return ((DBL)itd.Inv3 * (itd.pFn->Evaluate(VTmp) - threshold));
+    return ((DBL)itd.Inv3 * EvaluatePolarized (*itd.pFn, VTmp));
+}
+
+
+/*****************************************************************************/
+
+DBL IsoSurface::EvaluateAbs (GenericScalarFunctionInstance& fn, Vector3d& p) const
+{
+    return fabs (threshold - fn.Evaluate (p));
+}
+
+DBL IsoSurface::EvaluatePolarized (GenericScalarFunctionInstance& fn, Vector3d& p) const
+{
+    if (positivePolarity)
+        return threshold - fn.Evaluate (p);
+    else
+        return fn.Evaluate (p) - threshold;
+}
+
+bool IsoSurface::IsInside (GenericScalarFunctionInstance& fn, Vector3d& p) const
+{
+    if (positivePolarity)
+        return threshold < fn.Evaluate (p);
+    else
+        return fn.Evaluate (p) < threshold;
 }
 
 }