Commit 956b5591 authored by Will Schroeder's avatar Will Schroeder

Refactored kernel API to support different basis

parent eadf3ec1
Pipeline #7873 passed with stage
......@@ -3,6 +3,7 @@ set(Module_SRCS
vtkGaussianKernel.cxx
vtkInterpolationKernel.cxx
vtkPointInterpolator.cxx
vtkLinearKernel.cxx
vtkSPHKernel.cxx
vtkShepardKernel.cxx
vtkVoronoiKernel.cxx
......
......@@ -100,86 +100,84 @@ Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
//----------------------------------------------------------------------------
vtkIdType vtkEllipsoidalGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
ComputeBasis(double x[3], vtkIdList *pIds)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
vtkIdType numPts = pIds->GetNumberOfIds();
return pIds->GetNumberOfIds();
}
if ( numPts >= 1 ) //Use Gaussian kernel
//----------------------------------------------------------------------------
vtkIdType vtkEllipsoidalGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
vtkIdType numPts = pIds->GetNumberOfIds();
int i;
vtkIdType id;
double sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
double y[3], v[3], r2, z2, rxy2, mag;
double n[3], s;
double f2=this->F2, e2=this->E2;
for (i=0; i<numPts; ++i)
{
int i;
vtkIdType id;
double sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
double y[3], v[3], r2, z2, rxy2, mag;
double n[3], s;
double f2=this->F2, e2=this->E2;
for (i=0; i<numPts; ++i)
{
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
v[0] = x[0] - y[0];
v[1] = x[1] - y[1];
v[2] = x[2] - y[2];
r2 = vtkMath::Dot(v,v);
v[0] = x[0] - y[0];
v[1] = x[1] - y[1];
v[2] = x[2] - y[2];
r2 = vtkMath::Dot(v,v);
if ( r2 == 0.0 ) //precise hit on existing point
{
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
else // continue computing weights
{
// Normal affect
if ( this->Normals )
{
this->Normals->GetTuple(id,n);
mag = vtkMath::Dot(n,n);
mag = ( mag == 0.0 ? 1.0 : sqrt(mag) );
}
else
{
mag = 1.0;
}
if ( r2 == 0.0 ) //precise hit on existing point
// Scalar scaling
if ( this->Scalars )
{
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
this->Scalars->GetTuple(id,&s);
}
else // continue computing weights
else
{
// Normal affect
if ( this->Normals )
{
this->Normals->GetTuple(id,n);
mag = vtkMath::Dot(n,n);
mag = ( mag == 0.0 ? 1.0 : sqrt(mag) );
}
else
{
mag = 1.0;
}
// Scalar scaling
if ( this->Scalars )
{
this->Scalars->GetTuple(id,&s);
}
else
{
s = 1.0;
}
z2 = vtkMath::Dot(v,n) / mag;
z2 = z2*z2;
rxy2 = r2 - z2;
w[i] = s * exp(-f2 * (rxy2/e2 + z2));
sum += w[i];
}//computing weights
}//over all points
// Normalize
for (i=0; i<numPts; ++i)
{
w[i] /= sum;
}
s = 1.0;
}
return numPts;
}//using ellipsoidal Gaussian Kernel
z2 = vtkMath::Dot(v,n) / mag;
z2 = z2*z2;
rxy2 = r2 - z2;
else //null point
w[i] = s * exp(-f2 * (rxy2/e2 + z2));
sum += w[i];
}//computing weights
}//over all points
// Normalize
for (i=0; i<numPts; ++i)
{
return 0;
w[i] /= sum;
}
return numPts;
}
//----------------------------------------------------------------------------
......
......@@ -69,11 +69,21 @@ public:
vtkPointData *pd);
// Description:
// Given a point x, compute interpolation weights associated with nearby
// points. The method returns the number of nearby points N (i.e., the
// neighborhood). Note that both the nearby points list pIds and the
// weights array are of length N, are provided by the caller of the method,
// and may be dynamically resized as necessary.
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeWeights().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds);
// Description:
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are of length numPts, are
// provided by the caller of the method, and may be dynamically resized as
// necessary. Typically this method is called after ComputeBasis(),
// although advanced users can invoke ComputeWeights() and provide the
// interpolation basis points pIds directly.
virtual vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,
vtkDoubleArray *weights);
......
......@@ -50,54 +50,52 @@ Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
//----------------------------------------------------------------------------
vtkIdType vtkGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
ComputeBasis(double x[3], vtkIdList *pIds)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
vtkIdType numPts = pIds->GetNumberOfIds();
return pIds->GetNumberOfIds();
}
if ( numPts >= 1 ) //Use Gaussian kernel
//----------------------------------------------------------------------------
vtkIdType vtkGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
vtkIdType numPts = pIds->GetNumberOfIds();
int i;
vtkIdType id;
double d2, y[3], sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
double f2=this->F2;
for (i=0; i<numPts; ++i)
{
int i;
vtkIdType id;
double d2, y[3], sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
double f2=this->F2;
for (i=0; i<numPts; ++i)
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
d2 = vtkMath::Distance2BetweenPoints(x,y);
if ( d2 == 0.0 ) //precise hit on existing point
{
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
d2 = vtkMath::Distance2BetweenPoints(x,y);
if ( d2 == 0.0 ) //precise hit on existing point
{
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
else
{
w[i] = exp(-f2 * d2);
sum += w[i];
}
}//over all points
// Normalize
for (i=0; i<numPts; ++i)
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
else
{
w[i] /= sum;
w[i] = exp(-f2 * d2);
sum += w[i];
}
}//over all points
return numPts;
}//using Gaussian Kernel
else //null point
// Normalize
for (i=0; i<numPts; ++i)
{
return 0;
w[i] /= sum;
}
return numPts;
}
//----------------------------------------------------------------------------
......
......@@ -56,19 +56,29 @@ public:
// computational values.
virtual void Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds,
vtkPointData *pd);
// Description:
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeWeights().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds);
// Description:
// Given a point x, compute interpolation weights associated with nearby
// points. The method returns the number of nearby points N (i.e., the
// neighborhood). Note that both the nearby points list pIds and the
// weights array are of length N, are provided by the caller of the method,
// and may be dynamically resized as necessary.
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are provided by the caller
// of the method, and may be dynamically resized as necessary. Typically
// this method is called after ComputeBasis(), although advanced users can
// invoke ComputeWeights() and provide the interpolation basis points pIds
// directly.
virtual vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,
vtkDoubleArray *weights);
// Description:
// Specify the radius of the kernel. Points within this radius will be
// used for interpolation. If no point is found, then the closest point
// will be used.
// used for interpolation (if ComputeBasis() is invoked).
vtkSetClampMacro(Radius,double,0.000001,VTK_FLOAT_MAX);
vtkGetMacro(Radius,double);
......
......@@ -26,6 +26,13 @@
// method requires providing a point locator to accelerate neigborhood
// queries. Some kernels may refer back to the original dataset, or the point
// attribute data associated with the dataset.
//
// Typically the kernels are invoked in two parts: first, the basis is
// computed using the supplied point locator and dataset. These are the
// neighboring points used to compute the interpolation weights. Then, the
// weights are computed from points forming the basis. However, advanced
// users can develop their own basis, skipping the ComputeBasis() method, and
// then invoke ComputeWeights() directly.
// .SECTION Caveats
// The ComputeWeights() method is thread safe.
......@@ -77,11 +84,22 @@ public:
vtkBooleanMacro(RequiresInitialization, bool);
// Description:
// Given a point x, compute interpolation weights associated with nearby
// points. The method returns the number of nearby points N (i.e., the
// neighborhood). Note that both the nearby points list pIds and the
// weights array are of length N, are provided by the caller of the method,
// and may be dynamically resized as necessary.
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeBasis().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds) = 0;
// Description:
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are provided by the caller
// of the method, and may be dynamically resized as necessary. The method
// returns the number of weights (pIds may be resized in some
// cases). Typically this method is called after ComputeBasis(), although
// advanced users can invoke ComputeWeights() and provide the interpolation
// basis points pIds directly.
virtual vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,
vtkDoubleArray *weights) = 0;
......
/*=========================================================================
Program: Visualization Toolkit
Module: vtkLinearKernel.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 "vtkLinearKernel.h"
#include "vtkAbstractPointLocator.h"
#include "vtkObjectFactory.h"
#include "vtkIdList.h"
#include "vtkDoubleArray.h"
vtkStandardNewMacro(vtkLinearKernel);
//----------------------------------------------------------------------------
vtkLinearKernel::vtkLinearKernel()
{
this->Radius = 1.0;
}
//----------------------------------------------------------------------------
vtkLinearKernel::~vtkLinearKernel()
{
}
//----------------------------------------------------------------------------
vtkIdType vtkLinearKernel::
ComputeBasis(double x[3], vtkIdList *pIds)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
return pIds->GetNumberOfIds();
}
//----------------------------------------------------------------------------
vtkIdType vtkLinearKernel::
ComputeWeights(double*, vtkIdList *pIds, vtkDoubleArray *weights)
{
vtkIdType numPts = pIds->GetNumberOfIds();
double w = 1.0 / static_cast<double>(numPts);
weights->SetNumberOfTuples(numPts);
for (vtkIdType i=0; i < numPts; ++i)
{
weights->SetValue(i,w);
}
return numPts;
}
//----------------------------------------------------------------------------
void vtkLinearKernel::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkLinearKernel.h
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.
=========================================================================*/
// .NAME vtkLinearKernel - a linear interpolation kernel
// .SECTION Description
// vtkLinearKernel is an interpolation kernel that averages the contributions
// of all points in the basis.
// .SECTION See Also
// vtkInterpolationKernel vtkGaussianKernel vtkLinearKernel vtkShepardKernel
#ifndef vtkLinearKernel_h
#define vtkLinearKernel_h
#include "vtkFiltersPointsModule.h" // For export macro
#include "vtkInterpolationKernel.h"
class vtkIdList;
class vtkDoubleArray;
class VTKFILTERSPOINTS_EXPORT vtkLinearKernel : public vtkInterpolationKernel
{
public:
// Description:
// Standard methods for instantiation, obtaining type information, and printing.
static vtkLinearKernel *New();
vtkTypeMacro(vtkLinearKernel,vtkInterpolationKernel);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeWeights().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds);
// Description:
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are provided by the caller
// of the method, and may be dynamically resized as necessary. Typically
// this method is called after ComputeBasis(), although advanced users can
// invoke ComputeWeights() and provide the interpolation basis points pIds
// directly.
virtual vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,
vtkDoubleArray *weights);
// Description:
// Specify the radius of the kernel. Points within this radius will be used
// for interpolation (if ComputeBasis() is used).
vtkSetClampMacro(Radius,double,0.000001,VTK_FLOAT_MAX);
vtkGetMacro(Radius,double);
protected:
vtkLinearKernel();
~vtkLinearKernel();
double Radius;
private:
vtkLinearKernel(const vtkLinearKernel&); // Not implemented.
void operator=(const vtkLinearKernel&); // Not implemented.
};
#endif
......@@ -256,10 +256,9 @@ struct ProbePoints
{
this->Input->GetPoint(ptId,x);
numWeights = this->Kernel->ComputeWeights(x, pIds, weights);
if ( numWeights > 0)
if ( this->Kernel->ComputeBasis(x, pIds) > 0 )
{
numWeights = this->Kernel->ComputeWeights(x, pIds, weights);
this->Arrays.Interpolate(numWeights, pIds->GetPointer(0),
weights->GetPointer(0), ptId);
}
......@@ -324,10 +323,9 @@ struct ImageProbePoints : public ProbePoints
x[0] = origin[0] + i*spacing[0];
ptId = i + jOffset + kOffset;
numWeights = this->Kernel->ComputeWeights(x, pIds, weights);
if ( numWeights > 0)
if ( this->Kernel->ComputeBasis(x, pIds) > 0 )
{
numWeights = this->Kernel->ComputeWeights(x, pIds, weights);
this->Arrays.Interpolate(numWeights, pIds->GetPointer(0),
weights->GetPointer(0), ptId);
}
......
......@@ -34,11 +34,19 @@ vtkSPHKernel::~vtkSPHKernel()
//----------------------------------------------------------------------------
vtkIdType vtkSPHKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
ComputeBasis(double x[3], vtkIdList *pIds)
{
pIds->SetNumberOfIds(1);
vtkIdType pId = this->Locator->FindClosestPoint(x);
pIds->SetId(0,pId);
return 1;
}
//----------------------------------------------------------------------------
vtkIdType vtkSPHKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
......
......@@ -12,10 +12,10 @@
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
// .NAME vtkSPHKernel - a Voronoi interpolation kernel
// .NAME vtkSPHKernel - a SPH interpolation kernel
// .SECTION Description
// vtkSPHKernel is an interpolation kernel that simply returns the
// vtkSPHKernel is an interpolation kernel that uses a spherical h
// closest point to a point to be interpolated. A single weight ise returned
// with value=1.0.
......@@ -47,11 +47,21 @@ public:
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Given a point x, compute interpolation weights associated with nearby
// points. The method returns the number of nearby points N (i.e., the
// neighborhood). Note that both the nearby points list pIds and the
// weights array are of length N, are provided by the caller of the method,
// and may be dynamically resized as necessary.
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeWeights().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds);
// Description:
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are of length numPts, are
// provided by the caller of the method, and may be dynamically resized as
// necessary. Typically this method is called after ComputeBasis(),
// although advanced users can invoke ComputeWeights() and provide the
// interpolation basis points pIds directly.
virtual vtkIdType ComputeWeights(double x[3], vtkIdList *pIds,
vtkDoubleArray *weights);
......
......@@ -39,58 +39,57 @@ vtkShepardKernel::~vtkShepardKernel()
//----------------------------------------------------------------------------
vtkIdType vtkShepardKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
ComputeBasis(double x[3], vtkIdList *pIds)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
return pIds->GetNumberOfIds();
}
//----------------------------------------------------------------------------
vtkIdType vtkShepardKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
vtkIdType numPts = pIds->GetNumberOfIds();
int i;
vtkIdType id;
double d, y[3], sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
if ( numPts >= 1 ) //Use Shepard kernel
for (i=0; i<numPts; ++i)
{
int i;
vtkIdType id;
double d, y[3], sum = 0.0;
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
for (i=0; i<numPts; ++i)
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
if ( this->PowerParameter == 2.0 )
{
d = vtkMath::Distance2BetweenPoints(x,y);
}
else
{
id = pIds->GetId(i);
this->DataSet->GetPoint(id,y);
if ( this->PowerParameter == 2.0 )
{
d = vtkMath::Distance2BetweenPoints(x,y);
}
else
{
d = pow(sqrt(vtkMath::Distance2BetweenPoints(x,y)), this->PowerParameter);
}
if ( d == 0.0 ) //precise hit on existing point
{
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
else
{
w[i] = 1.0 / d;
sum += w[i];
}
}//over all points
// Normalize
for (i=0; i<numPts; ++i)
d = pow(sqrt(vtkMath::Distance2BetweenPoints(x,y)), this->PowerParameter);
}
if ( d == 0.0 ) //precise hit on existing point
{
w[i] /= sum;
pIds->SetNumberOfIds(1);
pIds->SetId(0,id);
weights->SetNumberOfTuples(1);
weights->SetValue(0,1.0);
return 1;
}
return numPts;
}//using Shepard Kernel
else
{
w[i] = 1.0 / d;
sum += w[i];
}
}//over all points
else //hit null point
// Normalize
for (i=0; i<numPts; ++i)
{
return 0;
w[i] /= sum;
}
return numPts;
}
//----------------------------------------------------------------------------
......
......@@ -50,11 +50,21 @@ public:
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Given a point x, compute interpolation weights associated with nearby
// points. The method returns the number of nearby points N (i.e., the
// neighborhood). Note that both the nearby points list pIds and the
// weights array are of length N, are provided by the caller of the method,
// and may be dynamically resized as necessary.
// Given a point x, determine the points around x which form an
// interpolation basis. The user must provide the vtkIdList pids, which will
// be dynamically resized as necessary. The method returns the number of
// points in the basis. Typically this method is called before
// ComputeWeights().
virtual vtkIdType ComputeBasis(double x[3], vtkIdList *pIds);
// Description:
// Given a point x, and a list of basis points pIds, compute interpolation
// weights associated with these basis points. Note that both the nearby
// basis points list pIds and the weights array are of length numPts, are
// provided by the caller of the method, and may be dynamically resized as