Commit c5e16cca authored by Bill Lorensen's avatar Bill Lorensen
Browse files

ENH: splines can now be closed (.i.e. periodic

parent 0e619d96
......@@ -41,9 +41,8 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
#include "vtkCardinalSpline.h"
// Description
// Construct a CardinalSpline wth the folloing defaults:
// ClampingOff
// Description:
// Construct a Cardinal Spline.
vtkCardinalSpline::vtkCardinalSpline ()
{
}
......@@ -66,6 +65,8 @@ float vtkCardinalSpline::Evaluate (float t)
intervals = this->Intervals;
coefficients = this->Coefficients;
if ( this->Closed ) size = size + 1;
// clamp the function at both ends
if (t < intervals[0]) t = intervals[0];
if (t > intervals[size - 1]) t = intervals[size - 1];
......@@ -101,48 +102,90 @@ void vtkCardinalSpline::Compute ()
// get the size of the independent variables
size = this->PiecewiseFunction->GetSize ();
// copy the independent variables
// copy the independent variables. Note that if the spline
// is closed the first and last point are assumed repeated -
// so we add and extra point
if (this->Intervals) delete [] this->Intervals;
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size; i++)
if ( !this->Closed )
{
this->Intervals[i] = *(ts + 2 * i);
}
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size; i++)
{
this->Intervals[i] = *(ts + 2 * i);
}
// allocate memory for work arrays
work = new float[size];
// allocate memory for work arrays
work = new float[size];
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for dependent variables
dependent = new float [size];
// allocate memory for dependent variables
dependent = new float [size];
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size; j++)
{
*(dependent + j) = *(xs + 2*j);
}
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size; j++)
{
*(dependent + j) = *(xs + 2*j);
}
this->Fit1D (size, this->Intervals, dependent,
this->Fit1D (size, this->Intervals, dependent,
work, (float (*)[4])coefficients,
this->LeftConstraint, this->LeftValue,
this->RightConstraint, this->RightValue);
}
else //add extra "fictitious" point to close loop
{
size = size + 1;
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size-1; i++)
{
this->Intervals[i] = *(ts + 2 * i);
}
this->Intervals[size-1] = this->Intervals[size-2] + 1.0;
// allocate memory for work arrays
work = new float[size];
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for dependent variables
dependent = new float [size];
// free the work array and dependent variable storage
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size-1; j++)
{
*(dependent + j) = *(xs + 2*j);
}
dependent[size-1] = *xs;
this->FitClosed1D (size, this->Intervals, dependent,
work, (float (*)[4])coefficients);
}
// free the work array and dependent variable storage
delete [] work;
delete [] dependent;
}
// Description
// Compute the coefficients for a 1D spline
// Compute the coefficients for a 1D spline. The spline is open.
void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
float *work, float coefficients[][4],
int leftConstraint, float leftValue,
......@@ -154,9 +197,8 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
int k;
// develop constraint at leftmost point.
switch (leftConstraint) {
switch (leftConstraint)
{
case 1:
// desired slope at leftmost point is leftValue.
coefficients[0][1] = 1.0;
......@@ -168,7 +210,7 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
coefficients[0][1] = 2.0;
coefficients[0][2] = 1.0;
work[0]= 3.0 * ((y[1] - y[0]) / (x[1] - x[0])) -
0.5 * (x[1] -x[0]) * leftValue;
0.5 * (x[1] - x[0]) * leftValue;
break;
case 3:
// desired second derivative at leftmost point is
......@@ -192,38 +234,38 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
coefficients[k][2] = xlk;
work[k] = 3.0 * (((xlkp * (y[k] - y[k-1])) / xlk) +
((xlk * (y[k+1] - y[k])) / xlkp));
}
}
// develop constraint at rightmost point.
switch (rightConstraint) {
case 1:
// desired slope at rightmost point is rightValue
coefficients[size - 1][0] = 0.0;
coefficients[size - 1][1] = 1.0;
work[size - 1] = rightValue;
break;
case 2:
// desired second derivative at rightmost point is rightValue.
coefficients[size-1][0] = 1.0;
coefficients[size-1][1] = 2.0;
work[size-1] = 3.0 * ((y[size-1] - y[size-2]) /
(x[size-1] - x[size-2])) +
0.5 * (x[size-1]-x[size-2]) * rightValue;
break;
case 3:
// desired second derivative at rightmost point is
// rightValue times second derivative at last interior point.
coefficients[size-1][0] = 4.0 * ((0.5 + rightValue) /
(2.0 + rightValue));
coefficients[size-1][1] = 2.0;
work[size-1] = 6.0 * ((1.0 + rightValue) / (2.0 + rightValue)) *
((y[size-1] - y[size-2]) /
(x[size-1] - x[size-2]));
break;
}
switch (rightConstraint)
{
case 1:
// desired slope at rightmost point is rightValue
coefficients[size - 1][0] = 0.0;
coefficients[size - 1][1] = 1.0;
work[size - 1] = rightValue;
break;
case 2:
// desired second derivative at rightmost point is rightValue.
coefficients[size-1][0] = 1.0;
coefficients[size-1][1] = 2.0;
work[size-1] = 3.0 * ((y[size-1] - y[size-2]) /
(x[size-1] - x[size-2])) +
0.5 * (x[size-1]-x[size-2]) * rightValue;
break;
case 3:
// desired second derivative at rightmost point is
// rightValue times second derivative at last interior point.
coefficients[size-1][0] = 4.0 * ((0.5 + rightValue) /
(2.0 + rightValue));
coefficients[size-1][1] = 2.0;
work[size-1] = 6.0 * ((1.0 + rightValue) / (2.0 + rightValue)) *
((y[size-1] - y[size-2]) /
(x[size-1] - x[size-2]));
break;
}
// solve resulting set of equations.
coefficients[0][2] = coefficients[0][2] / coefficients[0][1];
work[0] = work[0] / coefficients[0][1];
coefficients[size-1][2] = 0.0;
......@@ -246,7 +288,6 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
// derivative of the spline function at each joint.
// compute the coefficients of the cubic between
// each pair of joints.
for (k = 0; k < size - 1; k++)
{
b = x[k+1] - x[k];
......@@ -259,7 +300,7 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
}
// the coefficients of a fictitious nth cubic
// are evaluated. this may simplify
// are evaluated. This may simplify
// algorithms which include both end points.
coefficients[size-1][0] = y[size-1];
......@@ -270,6 +311,105 @@ void vtkCardinalSpline::Fit1D (int size, float *x, float *y,
}
// Description
// Compute the coefficients for a 1D spline. The spline is closed
// (i.e., the first and last point are assumed the same) and the
// spline is continous in value and derivatives.
void vtkCardinalSpline::FitClosed1D (int size, float *x, float *y,
float *work, float coefficients[][4])
{
float b;
float xlk;
float xlkp;
int k;
float aN, bN, cN, dN;
int N;
N = size - 1;
// develop body of band matrix.
//
for (k = 1; k < N; k++)
{
xlk = x[k] - x[k-1];
xlkp = x[k+1] - x[k];
coefficients[k][0] = xlkp;
coefficients[k][1] = 2.0 * (xlkp + xlk);
coefficients[k][2] = xlk;
work[k] = 3.0 * (((xlkp * (y[k] - y[k-1])) / xlk) +
((xlk * (y[k+1] - y[k])) / xlkp));
}
xlk = x[N] - x[N-1];
xlkp = x[1] - x[0];
aN = coefficients[N][0] = xlkp;
bN = coefficients[N][1] = 2.0 * (xlkp + xlk);
cN = coefficients[N][2] = xlk;
dN = work[N] = 3.0 * (((xlkp * (y[N] - y[N-1])) / xlk) +
((xlk * (y[1] - y[0])) / xlkp));
// solve resulting set of equations.
//
coefficients[0][2] = 0.0;
work[0] = 0.0;
coefficients[0][3] = 1.0;
for (k = 1; k <= N; k++)
{
coefficients[k][1] = coefficients[k][1] -
(coefficients[k][0] * coefficients[k-1][2]);
coefficients[k][2] = coefficients[k][2] / coefficients[k][1];
work[k] = (work[k] - (coefficients[k][0] * work[k-1]))
/ coefficients[k][1];
coefficients[k][3] = (-1.0 * coefficients[k][0] *
coefficients[k-1][3]) / coefficients[k][1];
}
coefficients[N][0] = 1.0;
coefficients[N][1] = 0.0;
for (k = N - 1; k > 0; k--)
{
coefficients[k][0] = coefficients[k][3] -
coefficients[k][2] * coefficients[k+1][0];
coefficients[k][1] = work[k] - coefficients[k][2] * coefficients[k+1][1];
}
work[0] = work[N] = (dN - cN * coefficients[1][1] -
aN * coefficients[N-1][1]) / ( bN +
cN * coefficients[1][0] + aN * coefficients[N-1][0]);
for (k=1; k < N; k++)
{
work[k] = coefficients[k][0] * work[N] + coefficients[k][1];
}
// the column vector work now contains the first
// derivative of the spline function at each joint.
// compute the coefficients of the cubic between
// each pair of joints.
for (k = 0; k < N; k++)
{
b = x[k+1] - x[k];
coefficients[k][0] = y[k];
coefficients[k][1] = work[k];
coefficients[k][2] = (3.0 * (y[k+1] - y[k])) / (b * b) -
(work[k+1] + 2.0 * work[k]) / b;
coefficients[k][3] = (2.0 * (y[k] - y[k+1])) / (b * b * b) +
(work[k+1] + work[k]) / (b * b);
}
// the coefficients of a fictitious nth cubic
// are the same as the coefficients in the first interval
//
coefficients[N][0] = y[N];
coefficients[N][1] = work[N];
coefficients[N][2] = coefficients[0][2];
coefficients[N][3] = coefficients[0][3];
}
void vtkCardinalSpline::PrintSelf(ostream& os, vtkIndent indent)
{
vtkSpline::PrintSelf(os,indent);
......
......@@ -41,14 +41,15 @@ MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
// .NAME vtkCardinalSpline - computes an interpolating spline using a
// a Cardinal basis.
// .SECTION Description
// vtkCardinalSpline is a concrete implementation of vtkSpline using a
// Cardinal basis.
// .SECTION See Also
// vtkSpline
// vtkSpline vtkKochanekSpline
#ifndef __vtkCardinalSpline_h
#define __vtkCardinalSpline_h
#include <stdio.h>
#include "vtkSpline.h"
class VTK_EXPORT vtkCardinalSpline : public vtkSpline
......@@ -68,6 +69,8 @@ public:
protected:
void Fit1D (int n, float *x, float *y, float *w, float coefficients[][4],
int leftConstraint, float leftValue, int rightConstraint, float rightValue);
void FitClosed1D (int n, float *x, float *y, float *w,
float coefficients[][4]);
};
#endif
......
......@@ -71,6 +71,8 @@ float vtkKochanekSpline::Evaluate (float t)
intervals = this->Intervals;
coefficients = this->Coefficients;
if ( this->Closed ) size = size + 1;
// clamp the function at both ends
if (t < intervals[0]) t = intervals[0];
if (t > intervals[size - 1]) t = intervals[size - 1];
......@@ -105,30 +107,64 @@ void vtkKochanekSpline::Compute ()
// get the size of the independent variables
size = this->PiecewiseFunction->GetSize ();
// copy the independent variables
if (this->Intervals) delete [] this->Intervals;
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size; i++)
if ( !this->Closed )
{
this->Intervals[i] = *(ts + 2 * i);
// copy the independent variables
if (this->Intervals) delete [] this->Intervals;
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size; i++)
{
this->Intervals[i] = *(ts + 2 * i);
}
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for dependent variables
dependent = new float [size];
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size; j++)
{
*(dependent + j) = *(xs + 2*j);
}
}
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for dependent variables
dependent = new float [size];
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size; j++)
else //spline is closed, create extra "fictitious" point
{
*(dependent + j) = *(xs + 2*j);
size = size + 1;
// copy the independent variables
if (this->Intervals) delete [] this->Intervals;
this->Intervals = new float[size];
ts = this->PiecewiseFunction->GetDataPointer ();
for (i = 0; i < size-1; i++)
{
this->Intervals[i] = *(ts + 2 * i);
}
this->Intervals[size-1] = this->Intervals[size-2] + 1.0;
// allocate memory for coefficients
if (this->Coefficients) delete [] this->Coefficients;
this->Coefficients = new float [4 * size];
// allocate memory for dependent variables
dependent = new float [size];
// get start of coefficients for this dependent variable
coefficients = this->Coefficients;
// get the dependent variable values
xs = this->PiecewiseFunction->GetDataPointer () + 1;
for (int j = 0; j < size-1; j++)
{
*(dependent + j) = *(xs + 2*j);
}
dependent[size-1] = *xs;
}
this->Fit1D (size, this->Intervals, dependent,
......@@ -166,15 +202,16 @@ void vtkKochanekSpline::Fit1D (int size, float *x, float *y,
if (size == 2)
{
// two points, set coeffieints for a straight line
coefficients[0][3] = 0.0;
coefficients[1][3] = 0.0;
coefficients[0][2] = 0.0;
coefficients[1][2] = 0.0;
coefficients[0][1] = (y[1] - y[0]) / (x[1] - x[0]);
coefficients[1][1] = coefficients[0][1];
coefficients[0][0] = y[0];
coefficients[1][0] = y[1];
// two points, set coefficients for a straight line
coefficients[0][3] = 0.0;
coefficients[1][3] = 0.0;
coefficients[0][2] = 0.0;
coefficients[1][2] = 0.0;
coefficients[0][1] = (y[1] - y[0]) / (x[1] - x[0]);
coefficients[1][1] = coefficients[0][1];
coefficients[0][0] = y[0];
coefficients[1][0] = y[1];
return;
}
N = size - 1;
......@@ -206,64 +243,90 @@ void vtkKochanekSpline::Fit1D (int size, float *x, float *y,
coefficients[0][0] = y[0];
coefficients[N][0] = y[N];
switch (leftConstraint) {
case 1:
// desired slope at leftmost point is leftValue
coefficients[0][1] = leftValue;
break;
case 2:
// desired second derivative at leftmost point is leftValue
coefficients[0][1] = (6*(y[1] - y[0]) - 2*coefficients[1][2] - leftValue)
/ 4.0;
break;
case 3:
// desired secord derivative at leftmost point is leftValue
// times secod derivative at first interior point
if ((leftValue > (-2.0 + EPSILON)) ||
(leftValue < (-2.0 - EPSILON)))
{
coefficients[0][1] = (3*(1 + leftValue)*(y[1] - y[0]) -
(1 + 2*leftValue)*coefficients[1][2])
/ (2 + leftValue);
}
else
{
coefficients[0][1] = 0.0;
}
break;
}
switch (rightConstraint)
if ( this->Closed ) //the curve is continuous and closed at P0=Pn
{
case 1:
// desired slope at rightmost point is rightValue
coefficients[N][2] = rightValue;
break;
case 2:
// desired second derivative at rightmost point is rightValue
coefficients[N][2] = (6*(y[N] - y[N-1]) - 2*coefficients[N-1][1] +
rightValue) / 4.0;
break;
case 3:
// desired secord derivative at rightmost point is rightValue
// times secord derivative at last interior point
if ((rightValue > (-2.0 + EPSILON)) ||
(rightValue < (-2.0 - EPSILON)))
{
coefficients[N][2] = (3*(1 + rightValue)*(y[N] - y[N-1]) -
(1 + 2*rightValue)*coefficients[N-1][1])
/ (2 + rightValue);
}
else
{
coefficients[N][2] = 0.0;
}
break;
cs = y[N] - y[N-1];
cd = y[1] - y[0];
ds = cs*((1 - tension)*(1 - continuity)*(1 + bias)) / 2.0
+ cd*((1 - tension)*(1 + continuity)*(1 - bias)) / 2.0;
dd = cs*((1 - tension)*(1 + continuity)*(1 + bias)) / 2.0
+ cd*((1 - tension)*(1 - continuity)*(1 - bias)) / 2.0;
// adjust deriviatives for non uniform spacing between nodes
n1 = x[1] - x[0];
n0 = x[N] - x[N-1];
ds *= (2 * n1 / (n0 + n1));
dd *= (2 * n0 / (n0 + n1));
coefficients[0][1] = dd;
coefficients[0][2] = ds;
coefficients[N][1] = dd;
coefficients[N][2] = ds;
}
else //curve is open
{
switch (leftConstraint)
{
case 1:
// desired slope at leftmost point is leftValue
coefficients[0][1] = leftValue;
break;
case 2:
// desired second derivative at leftmost point is leftValue
coefficients[0][1] = (6*(y[1] - y[0]) - 2*coefficients[1][2] - leftValue)
/ 4.0;
break;
case 3:
// desired secord derivative at leftmost point is leftValue
// times secod derivative at first interior point
if ((leftValue > (-2.0 + EPSILON)) ||
(leftValue < (-2.0 - EPSILON)))
{
coefficients[0][1] = (3*(1 + leftValue)*(y[1] - y[0]) -
(1 + 2*leftValue)*coefficients[1][2])
/ (2 + leftValue);