Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
953 lines (808 sloc) 23.4 KB
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPiecewiseFunction.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkPiecewiseFunction.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include <vtkstd/vector>
#include <vtkstd/set>
#include <vtkstd/algorithm>
#include <vtkstd/iterator>
vtkStandardNewMacro(vtkPiecewiseFunction);
// The Node structure
class vtkPiecewiseFunctionNode
{
public:
double X;
double Y;
double Sharpness;
double Midpoint;
};
// A comparison method for sorting nodes in increasing order
class vtkPiecewiseFunctionCompareNodes
{
public:
bool operator () ( const vtkPiecewiseFunctionNode *node1,
const vtkPiecewiseFunctionNode *node2 )
{
return node1->X < node2->X;
}
};
// A find method for finding a particular node in the function
class vtkPiecewiseFunctionFindNodeEqual
{
public:
double X;
bool operator () ( const vtkPiecewiseFunctionNode *node )
{
return node->X == this->X;
}
};
// A find method for finding nodes inside a specified range
class vtkPiecewiseFunctionFindNodeInRange
{
public:
double X1;
double X2;
bool operator () ( const vtkPiecewiseFunctionNode *node )
{
return ( node->X >= this->X1 &&
node->X <= this->X2 );
}
};
// A find method for finding nodes outside a specified range
class vtkPiecewiseFunctionFindNodeOutOfRange
{
public:
double X1;
double X2;
bool operator () ( const vtkPiecewiseFunctionNode *node )
{
return ( node->X < this->X1 ||
node->X > this->X2 );
}
};
// The internal structure for containing the STL objects
class vtkPiecewiseFunctionInternals
{
public:
vtkstd::vector<vtkPiecewiseFunctionNode*> Nodes;
vtkPiecewiseFunctionCompareNodes CompareNodes;
vtkPiecewiseFunctionFindNodeEqual FindNodeEqual;
vtkPiecewiseFunctionFindNodeInRange FindNodeInRange;
vtkPiecewiseFunctionFindNodeOutOfRange FindNodeOutOfRange;
};
// Construct a new vtkPiecewiseFunction with default values
vtkPiecewiseFunction::vtkPiecewiseFunction()
{
this->Clamping = 1;
this->Range[0] = 0;
this->Range[1] = 0;
this->Function = NULL;
this->AllowDuplicateScalars = 0;
this->Internal = new vtkPiecewiseFunctionInternals;
}
// Destruct a vtkPiecewiseFunction
vtkPiecewiseFunction::~vtkPiecewiseFunction()
{
if( this->Function )
{
delete [] this->Function;
}
for(unsigned int i=0;i<this->Internal->Nodes.size();i++)
{
delete this->Internal->Nodes[i];
}
this->Internal->Nodes.clear();
delete this->Internal;
}
void vtkPiecewiseFunction::DeepCopy( vtkDataObject *o )
{
vtkPiecewiseFunction *f = vtkPiecewiseFunction::SafeDownCast(o);
if (f != NULL)
{
this->Clamping = f->Clamping;
int i;
this->RemoveAllPoints();
for ( i = 0; i < f->GetSize(); i++ )
{
double val[4];
f->GetNodeValue(i, val);
this->AddPoint(val[0], val[1], val[2], val[3]);
}
this->Modified();
}
// Do the superclass
this->Superclass::DeepCopy(o);
}
void vtkPiecewiseFunction::ShallowCopy( vtkDataObject *o )
{
vtkPiecewiseFunction *f = vtkPiecewiseFunction::SafeDownCast(o);
if (f != NULL)
{
this->Clamping = f->Clamping;
int i;
this->RemoveAllPoints();
for ( i = 0; i < f->GetSize(); i++ )
{
double val[4];
f->GetNodeValue(i, val);
this->AddPoint(val[0], val[1], val[2], val[3]);
}
this->Modified();
}
// Do the superclass
this->vtkDataObject::ShallowCopy(o);
}
// This is a legacy method that is no longer needed
void vtkPiecewiseFunction::Initialize()
{
this->RemoveAllPoints();
}
// Return the number of points which specify this function
int vtkPiecewiseFunction::GetSize()
{
return static_cast<int>(this->Internal->Nodes.size());
}
// Return the type of function stored in object:
// Function Types:
// 0 : Constant (No change in slope between end points)
// 1 : NonDecreasing (Always increasing or zero slope)
// 2 : NonIncreasing (Always decreasing or zero slope)
// 3 : Varied (Contains both decreasing and increasing slopes)
// 4 : Unknown (Error condition)
//
const char *vtkPiecewiseFunction::GetType()
{
unsigned int i;
double value;
double prev_value = 0.0;
int function_type;
function_type = 0;
if( this->Internal->Nodes.size() )
{
prev_value = this->Internal->Nodes[0]->Y;
}
for( i=1; i < this->Internal->Nodes.size(); i++ )
{
value = this->Internal->Nodes[i]->Y;
// Do not change the function type if equal
if( value != prev_value )
{
if( value > prev_value )
{
switch( function_type )
{
case 0:
case 1:
function_type = 1; // NonDecreasing
break;
case 2:
function_type = 3; // Varied
break;
}
}
else // value < prev_value
{
switch( function_type )
{
case 0:
case 2:
function_type = 2; // NonIncreasing
break;
case 1:
function_type = 3; // Varied
break;
}
}
}
prev_value = value;
// Exit loop if we find a Varied function
if( function_type == 3 )
{
break;
}
}
switch( function_type )
{
case 0:
return "Constant";
case 1:
return "NonDecreasing";
case 2:
return "NonIncreasing";
case 3:
return "Varied";
}
return "Unknown";
}
// Since we no longer store the data in an array, we must
// copy out of the vector into an array. No modified check -
// could be added if performance is a problem
double *vtkPiecewiseFunction::GetDataPointer()
{
int size = static_cast<int>(this->Internal->Nodes.size());
if ( this->Function )
{
delete [] this->Function;
this->Function = NULL;
}
if ( size > 0 )
{
this->Function = new double[size*2];
for ( int i = 0; i < size; i++ )
{
this->Function[2*i ] = this->Internal->Nodes[i]->X;
this->Function[2*i+1] = this->Internal->Nodes[i]->Y;
}
}
return this->Function;
}
// Returns the first point location which starts a non-zero segment of the
// function. Note that the value at this point may be zero.
double vtkPiecewiseFunction::GetFirstNonZeroValue()
{
unsigned int i;
int all_zero = 1;
double x = 0.0;
// Check if no points specified
if( this->Internal->Nodes.size() == 0 )
{
return 0;
}
for( i=0; i < this->Internal->Nodes.size(); i++ )
{
if( this->Internal->Nodes[i]->Y != 0.0 )
{
x = this->Internal->Nodes[i]->X;
all_zero = 0;
break;
}
}
// If every specified point has a zero value then return
// a large value
if( all_zero )
{
x = VTK_DOUBLE_MAX;
}
else // A point was found with a non-zero value
{
if( i > 0 )
// Return the value of the point that precedes this one
{
x = this->Internal->Nodes[i-1]->X;
}
else
// If this is the first point in the function, return its
// value is clamping is off, otherwise VTK_DOUBLE_MIN if
// clamping is on.
{
if ( this->Clamping )
{
x = VTK_DOUBLE_MIN;
}
else
{
x = this->Internal->Nodes[0]->X;
}
}
}
return x;
}
// For a specified index value, get the node parameters
int vtkPiecewiseFunction::GetNodeValue( int index, double val[4] )
{
int size = static_cast<int>(this->Internal->Nodes.size());
if ( index < 0 || index >= size )
{
vtkErrorMacro("Index out of range!");
return -1;
}
val[0] = this->Internal->Nodes[index]->X;
val[1] = this->Internal->Nodes[index]->Y;
val[2] = this->Internal->Nodes[index]->Midpoint;
val[3] = this->Internal->Nodes[index]->Sharpness;
return 1;
}
// For a specified index value, get the node parameters
int vtkPiecewiseFunction::SetNodeValue( int index, double val[4] )
{
int size = static_cast<int>(this->Internal->Nodes.size());
if ( index < 0 || index >= size )
{
vtkErrorMacro("Index out of range!");
return -1;
}
this->Internal->Nodes[index]->X = val[0];
this->Internal->Nodes[index]->Y = val[1];
this->Internal->Nodes[index]->Midpoint = val[2];
this->Internal->Nodes[index]->Sharpness = val[3];
this->Modified();
return 1;
}
// Adds a point to the function. If a duplicate point is inserted
// then the function value at that location is set to the new value.
// This is the legacy version that assumes midpoint = 0.5 and
// sharpness = 0.0
int vtkPiecewiseFunction::AddPoint( double x, double y )
{
return this->AddPoint( x, y, 0.5, 0.0 );
}
// Adds a point to the function and returns the array index of the point.
int vtkPiecewiseFunction::AddPoint( double x, double y,
double midpoint, double sharpness )
{
// Error check
if ( midpoint < 0.0 || midpoint > 1.0 )
{
vtkErrorMacro("Midpoint outside range [0.0, 1.0]");
return -1;
}
if ( sharpness < 0.0 || sharpness > 1.0 )
{
vtkErrorMacro("Sharpness outside range [0.0, 1.0]");
return -1;
}
// remove any node already at this X location
if (!this->AllowDuplicateScalars)
{
this->RemovePoint( x );
}
// Create the new node
vtkPiecewiseFunctionNode *node = new vtkPiecewiseFunctionNode;
node->X = x;
node->Y = y;
node->Sharpness = sharpness;
node->Midpoint = midpoint;
// Add it, then sort to get everyting in order
this->Internal->Nodes.push_back(node);
this->SortAndUpdateRange();
// Now find this node so we can return the index
unsigned int i;
for ( i = 0; i < this->Internal->Nodes.size(); i++ )
{
if ( this->Internal->Nodes[i]->X == x )
{
break;
}
}
int retVal;
// If we didn't find it, something went horribly wrong so
// return -1
if ( i < this->Internal->Nodes.size() )
{
retVal = i;
}
else
{
retVal = -1;
}
return retVal;
}
// Sort the vector in increasing order, then fill in
// the Range
void vtkPiecewiseFunction::SortAndUpdateRange()
{
vtkstd::sort( this->Internal->Nodes.begin(),
this->Internal->Nodes.end(),
this->Internal->CompareNodes );
int size = static_cast<int>(this->Internal->Nodes.size());
if ( size )
{
this->Range[0] = this->Internal->Nodes[0]->X;
this->Range[1] = this->Internal->Nodes[size-1]->X;
}
else
{
this->Range[0] = 0;
this->Range[1] = 0;
}
this->Modified();
}
// Removes a point from the function. If no point is found then function
// remains the same.
int vtkPiecewiseFunction::RemovePoint( double x )
{
// First find the node since we need to know its
// index as our return value
unsigned int i;
for ( i = 0; i < this->Internal->Nodes.size(); i++ )
{
if ( this->Internal->Nodes[i]->X == x )
{
break;
}
}
int retVal;
// If the node doesn't exist, we return -1
if ( i < this->Internal->Nodes.size() )
{
retVal = i;
}
else
{
return -1;
}
// Now use STL to find it, so that we can remove it
this->Internal->FindNodeEqual.X = x;
vtkstd::vector<vtkPiecewiseFunctionNode*>::iterator iter =
vtkstd::find_if(this->Internal->Nodes.begin(),
this->Internal->Nodes.end(),
this->Internal->FindNodeEqual );
// Actually delete it
if ( iter != this->Internal->Nodes.end() )
{
delete *iter;
this->Internal->Nodes.erase(iter);
this->Modified();
}
else
{
// This should never happen - we already returned if the node
// didn't exist...
return -1;
}
return retVal;
}
// Removes all points from the function.
void vtkPiecewiseFunction::RemoveAllPoints()
{
for(unsigned int i=0;i<this->Internal->Nodes.size();i++)
{
delete this->Internal->Nodes[i];
}
this->Internal->Nodes.clear();
this->SortAndUpdateRange();
}
// Add in end points of line and remove any points between them
// Legacy method with no way to specify midpoint and sharpness
void vtkPiecewiseFunction::AddSegment( double x1, double y1,
double x2, double y2 )
{
int done;
// First, find all points in this range and remove them
done = 0;
while ( !done )
{
done = 1;
this->Internal->FindNodeInRange.X1 = x1;
this->Internal->FindNodeInRange.X2 = x2;
vtkstd::vector<vtkPiecewiseFunctionNode*>::iterator iter =
vtkstd::find_if(this->Internal->Nodes.begin(),
this->Internal->Nodes.end(),
this->Internal->FindNodeInRange );
if ( iter != this->Internal->Nodes.end() )
{
delete *iter;
this->Internal->Nodes.erase(iter);
this->Modified();
done = 0;
}
}
// Now add the points
this->AddPoint( x1, y1, 0.5, 0.0 );
this->AddPoint( x2, y2, 0.5, 0.0 );
}
// Return the value of the function at a position
double vtkPiecewiseFunction::GetValue( double x )
{
double table[1];
this->GetTable( x, x, 1, table );
return table[0];
}
// Remove all points outside the range, and make sure a point
// exists at each end of the range. Used as a convenience method
// for transfer function editors
int vtkPiecewiseFunction::AdjustRange(double range[2])
{
if (!range)
{
return 0;
}
double *function_range = this->GetRange();
// Make sure we have points at each end of the range
if (function_range[0] < range[0])
{
this->AddPoint(range[0], this->GetValue(range[0]));
}
else
{
this->AddPoint(range[0], this->GetValue(function_range[0]));
}
if (function_range[1] > range[1])
{
this->AddPoint(range[1], this->GetValue(range[1]));
}
else
{
this->AddPoint(range[1], this->GetValue(function_range[1]));
}
// Remove all points out-of-range
int done;
done = 0;
while ( !done )
{
done = 1;
this->Internal->FindNodeOutOfRange.X1 = range[0];
this->Internal->FindNodeOutOfRange.X2 = range[1];
vtkstd::vector<vtkPiecewiseFunctionNode*>::iterator iter =
vtkstd::find_if(this->Internal->Nodes.begin(),
this->Internal->Nodes.end(),
this->Internal->FindNodeOutOfRange );
if ( iter != this->Internal->Nodes.end() )
{
delete *iter;
this->Internal->Nodes.erase(iter);
this->Modified();
done = 0;
}
}
this->SortAndUpdateRange();
return 1;
}
// Returns a table of function values evaluated at regular intervals
void vtkPiecewiseFunction::GetTable( double xStart, double xEnd,
int size, double* table,
int stride )
{
int i;
int idx = 0;
int numNodes = static_cast<int>(this->Internal->Nodes.size());
// Need to keep track of the last value so that
// we can fill in table locations past this with
// this value if Clamping is On.
double lastValue = 0.0;
if ( numNodes != 0 )
{
lastValue = this->Internal->Nodes[numNodes-1]->Y;
}
double *tptr = NULL;
double x = 0.0;
double x1 = 0.0;
double x2 = 0.0;
double y1 = 0.0;
double y2 = 0.0;
double midpoint = 0.0;
double sharpness = 0.0;
// For each table entry
for ( i = 0; i < size; i++ )
{
// Find our location in the table
tptr = table + stride*i;
// Find our X location. If we are taking only 1 sample, make
// it halfway between start and end (usually start and end will
// be the same in this case)
if ( size > 1 )
{
x = xStart + (double(i)/double(size-1))*(xEnd-xStart);
}
else
{
x = 0.5*(xStart+xEnd);
}
// Do we need to move to the next node?
while ( idx < numNodes &&
x > this->Internal->Nodes[idx]->X )
{
idx++;
// If we are at a valid point index, fill in
// the value at this node, and the one before (the
// two that surround our current sample location)
// idx cannot be 0 since we just incremented it.
if ( idx < numNodes )
{
x1 = this->Internal->Nodes[idx-1]->X;
x2 = this->Internal->Nodes[idx ]->X;
y1 = this->Internal->Nodes[idx-1]->Y;
y2 = this->Internal->Nodes[idx ]->Y;
// We only need the previous midpoint and sharpness
// since these control this region
midpoint = this->Internal->Nodes[idx-1]->Midpoint;
sharpness = this->Internal->Nodes[idx-1]->Sharpness;
// Move midpoint away from extreme ends of range to avoid
// degenerate math
if ( midpoint < 0.00001 )
{
midpoint = 0.00001;
}
if ( midpoint > 0.99999 )
{
midpoint = 0.99999;
}
}
}
// Are we at the end? If so, just use the last value
if ( idx >= numNodes )
{
*tptr = (this->Clamping)?(lastValue):(0.0);
}
// Are we before the first node? If so, duplicate this nodes values
else if ( idx == 0 )
{
*tptr = (this->Clamping)?(this->Internal->Nodes[0]->Y):(0.0);
}
// Otherwise, we are between two nodes - interpolate
else
{
// Our first attempt at a normalized location [0,1] -
// we will be modifying this based on midpoint and
// sharpness to get the curve shape we want and to have
// it pass through (y1+y2)/2 at the midpoint.
double s = (x - x1) / (x2 - x1);
// Readjust based on the midpoint - linear adjustment
if ( s < midpoint )
{
s = 0.5 * s / midpoint;
}
else
{
s = 0.5 + 0.5*(s-midpoint)/(1.0-midpoint);
}
// override for sharpness > 0.99
// In this case we just want piecewise constant
if ( sharpness > 0.99 )
{
// Use the first value since we are below the midpoint
if ( s < 0.5 )
{
*tptr = y1;
continue;
}
// Use the second value at or above the midpoint
else
{
*tptr = y2;
continue;
}
}
// Override for sharpness < 0.01
// In this case we want piecewise linear
if ( sharpness < 0.01 )
{
// Simple linear interpolation
*tptr = (1-s)*y1 + s*y2;
continue;
}
// We have a sharpness between [0.01, 0.99] - we will
// used a modified hermite curve interpolation where we
// derive the slope based on the sharpness, and we compress
// the curve non-linearly based on the sharpness
// First, we will adjust our position based on sharpness in
// order to make the curve sharper (closer to piecewise constant)
if ( s < .5 )
{
s = 0.5 * pow(s*2,1.0 + 10*sharpness);
}
else if ( s > .5 )
{
s = 1.0 - 0.5 * pow((1.0-s)*2,1+10*sharpness);
}
// Compute some coefficients we will need for the hermite curve
double ss = s*s;
double sss = ss*s;
double h1 = 2*sss - 3*ss + 1;
double h2 = -2*sss + 3*ss;
double h3 = sss - 2*ss + s;
double h4 = sss - ss;
double slope;
double t;
// Use one slope for both end points
slope = y2 - y1;
t = (1.0 - sharpness)*slope;
// Compute the value
*tptr = h1*y1 + h2*y2 + h3*t + h4*t;
// Final error check to make sure we don't go outside
// the Y range
double min = (y1<y2)?(y1):(y2);
double max = (y1>y2)?(y1):(y2);
*tptr = (*tptr < min)?(min):(*tptr);
*tptr = (*tptr > max)?(max):(*tptr);
}
}
}
// Copy from double table to float
void vtkPiecewiseFunction::GetTable( double xStart, double xEnd,
int size, float* table,
int stride )
{
double *tmpTable = new double [size];
this->GetTable( xStart, xEnd, size, tmpTable, 1 );
double *tmpPtr = tmpTable;
float *tPtr = table;
for ( int i = 0; i < size; i++ )
{
*tPtr = static_cast<float>(*tmpPtr);
tPtr += stride;
tmpPtr ++;
}
delete[] tmpTable;
}
// Given a table of values, build the piecewise function. Legacy method
// that does not allow for midpoint and sharpness control
void vtkPiecewiseFunction::BuildFunctionFromTable( double xStart, double xEnd,
int size, double* table,
int stride )
{
double inc = 0.0;
double *tptr = table;
this->RemoveAllPoints();
if( size > 1 )
{
inc = (xEnd-xStart)/static_cast<double>(size-1);
}
int i;
for (i=0; i < size; i++)
{
vtkPiecewiseFunctionNode *node = new vtkPiecewiseFunctionNode;
node->X = xStart + inc*i;
node->Y = *tptr;
node->Sharpness = 0.0;
node->Midpoint = 0.5;
this->Internal->Nodes.push_back(node);
tptr += stride;
}
this->SortAndUpdateRange();
}
// Given a pointer to an array of values, build the piecewise function.
// Legacy method that does not allow for midpoint and sharpness control
void vtkPiecewiseFunction::FillFromDataPointer(int nb, double *ptr)
{
if (nb <= 0 || !ptr)
{
return;
}
this->RemoveAllPoints();
double *inPtr = ptr;
int i;
for (i=0; i < nb; i++)
{
vtkPiecewiseFunctionNode *node = new vtkPiecewiseFunctionNode;
node->X = inPtr[0];
node->Y = inPtr[1];
node->Sharpness = 0.0;
node->Midpoint = 0.5;
this->Internal->Nodes.push_back(node);
inPtr += 2;
}
this->SortAndUpdateRange();
}
//----------------------------------------------------------------------------
vtkPiecewiseFunction* vtkPiecewiseFunction::GetData(vtkInformation* info)
{
return
info? vtkPiecewiseFunction::SafeDownCast(info->Get(DATA_OBJECT())) : 0;
}
//----------------------------------------------------------------------------
vtkPiecewiseFunction* vtkPiecewiseFunction::GetData(vtkInformationVector* v,
int i)
{
return vtkPiecewiseFunction::GetData(v->GetInformationObject(i));
}
// Print method for tkPiecewiseFunction
void vtkPiecewiseFunction::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
unsigned int i;
os << indent << "Clamping: " << this->Clamping << endl;
os << indent << "Range: [" << this->Range[0] << ","
<< this->Range[1] << "]" << endl;
os << indent << "Function Points: " << this->Internal->Nodes.size() << endl;
for( i = 0; i < this->Internal->Nodes.size(); i++ )
{
os << indent << " " << i << " X: "
<< this->Internal->Nodes[i]->X << " Y: "
<< this->Internal->Nodes[i]->Y << " Sharpness: "
<< this->Internal->Nodes[i]->Sharpness << " Midpoint: "
<< this->Internal->Nodes[i]->Midpoint << endl;
}
os << indent << "AllowDuplicateScalars: " << this->AllowDuplicateScalars
<< endl;
}