Commit e5c4dd6d authored by Will Schroeder's avatar Will Schroeder

Added new class for oriented gaussian distributions

parent b92e2cbc
Pipeline #7876 passed with stage
......@@ -22,6 +22,7 @@ set(Module_SRCS
vtkDelaunay2D.cxx
vtkDelaunay3D.cxx
vtkElevationFilter.cxx
vtkEllipsoidalGaussianKernel.cxx
vtkExecutionTimer.cxx
vtkFeatureEdges.cxx
vtkFieldDataToAttributeDataFilter.cxx
......
......@@ -69,6 +69,7 @@ outlineActor.SetMapper(outlineMapper)
# Gaussian kernel-------------------------------------------------------
gaussianKernel = vtk.vtkGaussianKernel()
#gaussianKernel = vtk.vtkEllipsoidalGaussianKernel()
gaussianKernel.SetSharpness(4);
interpolator1 = vtk.vtkPointInterpolator()
......
/*=========================================================================
Program: Visualization Toolkit
Module: vtkEllipsoidalGaussianKernel.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 "vtkEllipsoidalGaussianKernel.h"
#include "vtkAbstractPointLocator.h"
#include "vtkObjectFactory.h"
#include "vtkIdList.h"
#include "vtkDoubleArray.h"
#include "vtkDataSet.h"
#include "vtkPointData.h"
#include "vtkMath.h"
vtkStandardNewMacro(vtkEllipsoidalGaussianKernel);
//----------------------------------------------------------------------------
vtkEllipsoidalGaussianKernel::vtkEllipsoidalGaussianKernel()
{
this->Radius = 1.0;
this->Sharpness = 2.0;
this->Eccentricity = 2.0;
this->Normals = NULL;
this->Scalars = NULL;
this->F2 = this->Sharpness / this->Radius;
this->E2 = this->Eccentricity * this->Eccentricity;
}
//----------------------------------------------------------------------------
vtkEllipsoidalGaussianKernel::~vtkEllipsoidalGaussianKernel()
{
this->FreeStructures();
}
//----------------------------------------------------------------------------
void vtkEllipsoidalGaussianKernel::
FreeStructures()
{
this->Superclass::FreeStructures();
if ( this->Normals )
{
this->Normals->Delete();
this->Normals = NULL;
}
if ( this->Scalars )
{
this->Scalars->Delete();
this->Scalars = NULL;
}
}
//----------------------------------------------------------------------------
void vtkEllipsoidalGaussianKernel::
Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
{
this->Superclass::Initialize(loc, ds, pd);
this->Scalars = pd->GetScalars();
if ( this->Scalars && this->Scalars->GetNumberOfComponents() == 1 )
{
this->Scalars->Register(this);
}
else
{
this->Scalars = NULL;
}
this->Normals = pd->GetNormals();
if ( this->Normals )
{
this->Normals->Register(this);
}
else
{
this->Normals = NULL;
}
this->F2 = this->Sharpness / this->Radius;
this->F2 = this->F2 * this->F2;
this->E2 = this->Eccentricity * this->Eccentricity;
}
//----------------------------------------------------------------------------
vtkIdType vtkEllipsoidalGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
vtkIdType numPts = pIds->GetNumberOfIds();
if ( numPts >= 1 ) //Use Gaussian kernel
{
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);
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;
}
// 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;
}
return numPts;
}//using ellipsoidal Gaussian Kernel
else //null point
{
return 0;
}
}
//----------------------------------------------------------------------------
void vtkEllipsoidalGaussianKernel::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Radius: " << this->Radius << endl;
os << indent << "Sharpness: " << this->Sharpness << endl;
os << indent << "Eccentricity: " << this->Eccentricity << endl;
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkEllipsoidalGaussianKernel.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 vtkEllipsoidalGaussianKernel - an ellipsoidal Gaussian interpolation kernel
// .SECTION Description
// vtkEllipsoidalGaussianKernel is an interpolation kernel that returns the
// weights for all points found in the ellipsoid defined by radius R in
// combination with local data (normals and/or scalars). For example, "pancake"
// weightings (the local normal parallel to the minimum ellisoidal axis); or
// "needle" weightings (the local normal parallel to the maximum ellipsoidal
// axis) are possible. (Note that spherical Gaussian weightings are more
// efficiently computed using vtkGaussianKernel.)
//
// The ellipsoidal Gaussian can be described by:
//
// W(x) = S * exp( -( Sharpness^2 * ((rxy/E)**2 + z**2)/R**2) )
//
// where S is the local scalar value; E is a user-defined eccentricity factor
// that controls the elliptical shape of the splat; z is the distance of the
// current voxel sample point along the local normal N; and rxy is the
// distance to neigbor point x in the direction prependicular to N.
// .SECTION Caveats
// The weights are normalized so that SUM(Wi) = 1. If a neighbor point p
// precisely lies on the point to be interpolated, then the interpolated
// point takes on the values associated with p.
// .SECTION See Also
// vtkPointInterpolator vtkInterpolationKernel vtkGaussianKernel
// vtkVoronoiKernel vtkSPHKernel vtkShepardKernel
#ifndef vtkEllipsoidalGaussianKernel_h
#define vtkEllipsoidalGaussianKernel_h
#include "vtkFiltersCoreModule.h" // For export macro
#include "vtkInterpolationKernel.h"
class vtkIdList;
class vtkDataArray;
class vtkDoubleArray;
class VTKFILTERSCORE_EXPORT vtkEllipsoidalGaussianKernel : public vtkInterpolationKernel
{
public:
// Description:
// Standard methods for instantiation, obtaining type information, and printing.
static vtkEllipsoidalGaussianKernel *New();
vtkTypeMacro(vtkEllipsoidalGaussianKernel,vtkInterpolationKernel);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Initialize the kernel. Overload the superclass to set up scalars and
// vectors.
virtual void Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds,
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.
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.
vtkSetClampMacro(Radius,double,0.000001,VTK_FLOAT_MAX);
vtkGetMacro(Radius,double);
// Description:
// Set / Get the sharpness (i.e., falloff) of the Gaussian. By default
// Sharpness=2. As the sharpness increases the effects of distant points
// are reduced.
vtkSetClampMacro(Sharpness,double,1,VTK_FLOAT_MAX);
vtkGetMacro(Sharpness,double);
// Description:
// Set / Get the eccentricity of the ellipsoidal Gaussian. A value=1.0
// produces a spherical distribution. Values < 1 produce a needle like
// distribution (in the direction of the normal); values > 1 produce a
// pancake like distribution (orthogonal to the normal).
vtkSetClampMacro(Eccentricity,double,0.000001,VTK_FLOAT_MAX);
vtkGetMacro(Eccentricity,double);
protected:
vtkEllipsoidalGaussianKernel();
~vtkEllipsoidalGaussianKernel();
double Radius;
double Sharpness;
double Eccentricity;
vtkDataArray *Normals;
vtkDataArray *Scalars;
// Internal structure to reduce computation
double F2, E2;
virtual void FreeStructures();
private:
vtkEllipsoidalGaussianKernel(const vtkEllipsoidalGaussianKernel&); // Not implemented.
void operator=(const vtkEllipsoidalGaussianKernel&); // Not implemented.
};
#endif
......@@ -28,6 +28,8 @@ vtkGaussianKernel::vtkGaussianKernel()
{
this->Radius = 1.0;
this->Sharpness = 2.0;
this->F2 = this->Sharpness / this->Radius;
}
......@@ -36,6 +38,15 @@ vtkGaussianKernel::~vtkGaussianKernel()
{
}
//----------------------------------------------------------------------------
void vtkGaussianKernel::
Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
{
this->Superclass::Initialize(loc, ds, pd);
this->F2 = this->Sharpness / this->Radius;
this->F2 = this->F2 * this->F2;
}
//----------------------------------------------------------------------------
vtkIdType vtkGaussianKernel::
......@@ -49,9 +60,9 @@ ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
int i;
vtkIdType id;
double d2, y[3], sum = 0.0;
double h2 = (this->Radius*this->Radius) / (this->Sharpness*this->Sharpness);
weights->SetNumberOfTuples(numPts);
double *w = weights->GetPointer(0);
double f2=this->F2;
for (i=0; i<numPts; ++i)
{
......@@ -69,7 +80,7 @@ ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
}
else
{
w[i] = exp(-(d2/h2));
w[i] = exp(-f2 * d2);
sum += w[i];
}
}//over all points
......@@ -94,4 +105,6 @@ void vtkGaussianKernel::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Radius: " << this->Radius << endl;
os << indent << "Sharpness: " << this->Sharpness << endl;
}
......@@ -16,10 +16,11 @@
// .SECTION Description
// vtkGaussianKernel is an interpolation kernel that simply returns the
// weights for all points found in the sphere defined by radius R. The weights
// are computed as: exp(-(s*r/R)^2) where r is the distance from the point to be
// interpolated to a neighboring point within R. The sharpness s simply affects
// the rate of fall off of the Gaussian.
// weights for all points found in the sphere defined by radius R. The
// weights are computed as: exp(-(s*r/R)^2) where r is the distance from the
// point to be interpolated to a neighboring point within R. The sharpness s
// simply affects the rate of fall off of the Gaussian. (A more general
// Gaussian kernel is available from vtkEllipsoidalGaussianKernel.)
// .SECTION Caveats
// The weights are normalized sp that SUM(Wi) = 1. If a neighbor point p
......@@ -27,8 +28,8 @@
// point takes on the values associated with p.
// .SECTION See Also
// vtkPointInterpolator vtkInterpolationKernel vtkVoronoiKernel vtkSPHKernel
// vtkShepardKernel
// vtkPointInterpolator vtkInterpolationKernel vtkEllipsoidalGaussianKernel
// vtkVoronoiKernel vtkSPHKernel vtkShepardKernel
#ifndef vtkGaussianKernel_h
......@@ -50,6 +51,11 @@ public:
vtkTypeMacro(vtkGaussianKernel,vtkInterpolationKernel);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Initialize the kernel. Overload the superclass to set up internal
// computational values.
virtual void Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds,
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
......@@ -63,7 +69,7 @@ public:
// 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.
vtkSetClampMacro(Radius,double,0.0,VTK_FLOAT_MAX);
vtkSetClampMacro(Radius,double,0.000001,VTK_FLOAT_MAX);
vtkGetMacro(Radius,double);
// Description:
......@@ -80,6 +86,9 @@ protected:
double Radius;
double Sharpness;
// Internal structure to reduce computation
double F2;
private:
vtkGaussianKernel(const vtkGaussianKernel&); // Not implemented.
void operator=(const vtkGaussianKernel&); // Not implemented.
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment