This function returns an estimate of the gradient of the entropy of a particle distribution with respect to change in position of a specific particle in that distribution. More...

#include <itkPSMParticleEntropyFunction.h>

Inheritance diagram for itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >:

Collaboration diagram for itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >:

Public Types
typedef PSMParticleEntropyFunction	Self

typedef SmartPointer< Self >	Pointer

typedef SmartPointer< const Self >	ConstPointer

typedef PSMCostFunction< VDimension >	Superclass

typedef TGradientNumericType	GradientNumericType

typedef Superclass::ParticleSystemType	ParticleSystemType

typedef PSMContainerArrayAttribute< double, VDimension >	SigmaCacheType

typedef Superclass::VectorType	VectorType

typedef ParticleSystemType::PointType	PointType

typedef vnl_vector_fixed< TGradientNumericType, VDimension >	GradientVectorType

Public Types inherited from itk::PSMCostFunction< VDimension >
typedef PSMCostFunction	Self

typedef SmartPointer< Self >	Pointer

typedef SmartPointer< const Self >	ConstPointer

typedef LightObject	Superclass

typedef PSMParticleSystem< VDimension >	ParticleSystemType

typedef vnl_vector_fixed< double, VDimension >	VectorType

Public Member Functions
	itkTypeMacro (PSMParticleEntropyFunction, PSMCostFunction)

	itkNewMacro (Self)

	itkStaticConstMacro (Dimension, unsigned int, VDimension)

virtual VectorType	Evaluate (unsigned int, unsigned int, const ParticleSystemType *, double &) const

virtual VectorType	Evaluate (unsigned int, unsigned int, const ParticleSystemType *, double &, double &) const

virtual double	Energy (unsigned int, unsigned int, const ParticleSystemType *) const

virtual void	ResetBuffers ()

virtual double	EstimateSigma (unsigned int, const typename ParticleSystemType::PointVectorType &, const std::vector< double > &, const PointType &, double, double, int &err) const

TGradientNumericType	AngleCoefficient (const GradientVectorType &, const GradientVectorType &) const

void	SetMinimumNeighborhoodRadius (double s)

double	GetMinimumNeighborhoodRadius () const

void	SetMaximumNeighborhoodRadius (double s)

double	GetMaximumNeighborhoodRadius () const

void	SetFlatCutoff (double s)

double	GetFlatCutoff () const

void	SetNeighborhoodToSigmaRatio (double s)

double	GetNeighborhoodToSigmaRatio () const

void	SetSpatialSigmaCache (SigmaCacheType *s)

SigmaCacheType *	GetSpatialSigmaCache ()

const SigmaCacheType *	GetSpatialSigmaCache () const

void	ComputeAngularWeights (const PointType &, const typename ParticleSystemType::PointVectorType &, const PSMImageDomainWithGradients< TGradientNumericType, VDimension > *, std::vector< double > &) const

Public Member Functions inherited from itk::PSMCostFunction< VDimension >
	itkTypeMacro (PSMCostFunction, LightObject)

	itkStaticConstMacro (Dimension, unsigned int, VDimension)

virtual void	AfterIteration ()

virtual void	BeforeIteration ()

virtual void	BeforeEvaluate (unsigned int, unsigned int, const ParticleSystemType *)

virtual void	SetParticleSystem (ParticleSystemType *p)

virtual ParticleSystemType *	GetParticleSystem () const

virtual void	SetDomainNumber (unsigned int i)

virtual int	GetDomainNumber () const

Protected Member Functions
void	operator= (const PSMParticleEntropyFunction &)

	PSMParticleEntropyFunction (const PSMParticleEntropyFunction &)

Protected Member Functions inherited from itk::PSMCostFunction< VDimension >
void	operator= (const PSMCostFunction &)

	PSMCostFunction (const PSMCostFunction &)

Protected Attributes
double	m_MinimumNeighborhoodRadius

double	m_MaximumNeighborhoodRadius

double	m_FlatCutoff

double	m_NeighborhoodToSigmaRatio

SigmaCacheType::Pointer	m_SpatialSigmaCache

Protected Attributes inherited from itk::PSMCostFunction< VDimension >
ParticleSystemType *	m_ParticleSystem

unsigned int	m_DomainNumber

Detailed Description

template<class TGradientNumericType, unsigned int VDimension>
class itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >

This function returns an estimate of the gradient of the entropy of a particle distribution with respect to change in position of a specific particle in that distribution.

The following description is an excerpt from

J Cates, P T Fletcher, M Styner, M Shenton, R Whitaker. Shape Modeling and Analysis with Entropy-Based Particle Systems. Information Processing in Medical Imaging IPMI 2007, LNCS 4584, pp. 333�345, 2007.

We treat a surface as a subset of $\Re^d$ , where $d=2$ or $d=3$ depending whether we are processing curves in the plane or surfaces in a volume, refspectively. The method we describe here deals with smooth, closed manifolds of codimension one, and we will refer to such manifolds as {surfaces}. We sample a surface ${\cal S} \subset \Re^d$ using a discrete set of $N$ points that are considered random variables $Z = (X_1, X_2, \ldots, X_N)$ drawn from a probability density function (PDF), $p(X)$ . We denote a realization of this PDF with lower case, and thus we have $z = (x_1, x_2,\ldots, x_N)$ , where $z \in {\cal S}^N$ . The probability of a realization $x$ is $p(X = x)$ , which we denote simply as $p(x)$ .

The amount of information contained in such a random sampling is, in the limit, the differential entropy of the PDF, which is $H[X] = -\int_S p(x) \log p(x) dx = -E\{\log p(X)\}$ , where $E\{ \cdot \}$ is the expectation. When we have a sufficient number of points sampled from $p$ , we can approximate the expectation by the sample mean, which gives $H[X] \approx - (1/N)\sum_{i} \log p(x_i)$ . We must also estimate $p(x_i)$ . Density functions on surfaces can be quite complex, and so we use a nonparametric, Parzen windowing estimation of this density using the particles themselves. Thus we have

$p(x_i) \approx \frac{1}{N(N-1)} \sum^N_{j=1, j \neq i} G(x_i - x_j, \sigma_i),$

where $G(x_i - x_j, \sigma_i)$ is a $d$ -dimensional, isotropic Gaussian with standard deviation $\sigma_i$ . The cost function $C$ , is therefore an approximation of (negative) entropy:

$-H[X] \approx C(x_1, \dots, x_N) = \sum_{i} \log \frac{1}{N(N-1)} \sum_{j \neq i} G(x_i - x_j, \sigma_i).$

Author: Josh Cates

Definition at line 84 of file itkPSMParticleEntropyFunction.h.

Member Typedef Documentation

template<class TGradientNumericType, unsigned int VDimension>

typedef TGradientNumericType itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::GradientNumericType

Data type representing individual gradient components.

Definition at line 95 of file itkPSMParticleEntropyFunction.h.

template<class TGradientNumericType, unsigned int VDimension>

typedef Superclass::ParticleSystemType itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::ParticleSystemType

Type of particle system.

Definition at line 98 of file itkPSMParticleEntropyFunction.h.

template<class TGradientNumericType, unsigned int VDimension>

typedef PSMParticleEntropyFunction itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::Self

Standard class typedefs.

Definition at line 88 of file itkPSMParticleEntropyFunction.h.

template<class TGradientNumericType, unsigned int VDimension>

typedef PSMContainerArrayAttribute<double, VDimension> itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SigmaCacheType

Cache type for the sigma values.

Definition at line 101 of file itkPSMParticleEntropyFunction.h.

template<class TGradientNumericType, unsigned int VDimension>

typedef Superclass::VectorType itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::VectorType

Vector & Point types.

Definition at line 104 of file itkPSMParticleEntropyFunction.h.

Member Function Documentation

template<class TGradientNumericType , unsigned int VDimension>

TGradientNumericType itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::AngleCoefficient	(	const GradientVectorType &	p_i_normal,
		const GradientVectorType &	p_j_normal
	)		const

Returns a weighting coefficient based on the angle between two vectors. Weights smoothly approach zero as the angle between two normals approaches 90 degrees.

Definition at line 27 of file itkPSMParticleEntropyFunction.hxx.

 {
   // Get the cosine of the angle between the two particles' normals
   TGradientNumericType cosine = dot_product(p_i_normal,p_j_normal) /
     (p_i_normal.magnitude()*p_j_normal.magnitude() + 1.0e-6);
   
   // the flat region
   if ( cosine >= m_FlatCutoff ) return 1.0;
   
   // the feathered region
   return ( cos( (m_FlatCutoff - cosine) / (1.0+m_FlatCutoff) * (3.14159265358979/2.0) )) ; 
 } 

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::ComputeAngularWeights	(	const PointType &	pos,
		const typename ParticleSystemType::PointVectorType &	neighborhood,
		const PSMImageDomainWithGradients< TGradientNumericType, VDimension > *	domain,
		std::vector< double > &	weights
	)		const

Compute a set of weights based on the difference in the normals of a central point and each of its neighbors. Difference of > 90 degrees results in a weight of 0.

Definition at line 43 of file itkPSMParticleEntropyFunction.hxx.

 {
   GradientVectorType posnormal = domain->SampleNormalVnl(pos, 1.0e-10);
   weights.resize(neighborhood.size());
   
   for (unsigned int i = 0; i < neighborhood.size(); i++)
     {
     weights[i] = this->AngleCoefficient(posnormal,
                                         domain->SampleNormalVnl(neighborhood[i].Point, 1.0e-10));
     if (weights[i] < 1.0e-5) weights[i] = 0.0;
     }
 }

template<class TGradientNumericType , unsigned int VDimension>

double itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::EstimateSigma	(	unsigned	int,
		const typename ParticleSystemType::PointVectorType &	neighborhood,
		const std::vector< double > &	weights,
		const PointType &	pos,
		double	initial_sigma,
		double	precision,
		int &	err
	)		const

virtual

Estimate the best sigma for Parzen windowing in a given neighborhood. The best sigma is the sigma that maximizes probability at the given point

Definition at line 62 of file itkPSMParticleEntropyFunction.hxx.

 {
   const double epsilon = 1.0e-5;
   const double min_sigma = 1.0e-4;
   
   const double M = static_cast<double>(VDimension);
   const double MM = M * M * 2.0 + M;
   
   double error = 1.0e6;
   double sigma, prev_sigma;
   sigma = initial_sigma;
   
   while (error > precision)
     {
     VectorType r_vec;
     double A = 0.0;
     double B = 0.0;
     double C = 0.0;
     double sigma2 = sigma * sigma;
     double sigma22 = sigma2 * 2.0;
     
     for (unsigned int i = 0; i < neighborhood.size(); i++)
       {
       if (weights[i] < epsilon) continue;
       
       //    if ( neighborhood[i].Index == idx) continue;
       for (unsigned int n = 0; n < VDimension; n++)
         {
         r_vec[n] = pos[n] - neighborhood[i].Point[n];
         }
       
       double r = r_vec.magnitude();
       double r2 = r*r;
       double alpha = exp(-r2 / sigma22) * weights[i];
       A += alpha;
       B += r2 * alpha;
       C += r2 * r2 * alpha;
       } // end for i
 
     prev_sigma = sigma;
     
     if (A < epsilon)
       {
       err = 1;
       return sigma;
       }; // results are not meaningful
     
     // First order convergence update.  This is a fixed point iteration.
     //sigma = sqrt(( 1.0 / DIM ) * ( B / A));
     
     // Second order convergence update (Newton-Raphson).  This is the first
     // derivative of the negative of the probability density estimation function squared over the
     // second derivative.
     
     // old math
     //    sigma -= (sigma2 * VDimension * A * A - A  * B) / (((2.0 * sigma * VDimension) * A * A -
     //                                          (1.0/(sigma2*sigma))*(A*C-B*B)) + epsilon);
     
     // New math -- results are not obviously different? 
     sigma -= (A * (B - A * sigma2 * M)) /
       ( (-MM * A *A * sigma) - 3.0 * A * B * (1.0 / (sigma + epsilon))
         - (A*C + B*B) * (1.0 / (sigma2 * sigma + epsilon)) + epsilon);
     
     error = 1.0 - fabs((sigma/prev_sigma));
     
     // Constrain sigma.
     if (sigma < min_sigma)
       {
       sigma = min_sigma;
       error = precision; // we are done if sigma has vanished
       }
     else
       {
       if (sigma < 0.0) sigma = min_sigma;
       }
     
     } // end while (error > precision)
   
   err = 0;
   return sigma;
   
 } // end estimate sigma

template<class TGradientNumericType , unsigned int VDimension>

PSMParticleEntropyFunction< TGradientNumericType, VDimension >::VectorType itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::Evaluate	(	unsigned int	idx,
		unsigned int	d,
		const ParticleSystemType *	system,
		double &	maxdt
	)		const

virtual

The first argument is a pointer to the particle system. The second argument is the index of the domain within that particle system. The third argument is the index of the particle location within the given domain.

Implements itk::PSMCostFunction< VDimension >.

Definition at line 152 of file itkPSMParticleEntropyFunction.hxx.

 {
   // Grab a pointer to the domain.  We need a Domain that has surface normal information.
   const PSMImageDomainWithGradients<TGradientNumericType, VDimension> *
     domain = static_cast<const PSMImageDomainWithGradients<
   TGradientNumericType, VDimension> *>(system->GetDomain(d));
   const double epsilon = 1.0e-6;
   
   // Retrieve the previous optimal sigma value for this point.  If the value is
   // tiny (i.e. unitialized) then use a fraction of the maximum allowed
   // neighborhood radius.
   double sigma = m_SpatialSigmaCache->operator[](d)->operator[](idx);
   if (sigma < epsilon)
     { sigma = m_MinimumNeighborhoodRadius / m_NeighborhoodToSigmaRatio;}
   
   // Determine the extent of the neighborhood that will be used in the Parzen
   // windowing estimation.  The neighborhood extent is based on the optimal
   // sigma calculation and limited to a user supplied maximum radius (probably
   // the size of the domain).
   double neighborhood_radius = sigma * m_NeighborhoodToSigmaRatio;
   if (neighborhood_radius > m_MaximumNeighborhoodRadius)
     { neighborhood_radius = m_MaximumNeighborhoodRadius; }
   
   // Get the position for which we are computing the gradient.
   PointType pos = system->GetPosition(idx, d);
   
   // Get the neighborhood surrounding the point "pos".
   typename ParticleSystemType::PointVectorType neighborhood
     = system->FindNeighborhoodPoints(pos, neighborhood_radius, d);
   
   // Compute the weights based on angle between the neighbors and the center.
   std::vector<double> weights;
   this->ComputeAngularWeights(pos,neighborhood,domain,weights);
   
   // Estimate the best sigma for Parzen windowing.  In some cases, such as when
   // the neighborhood does not include enough points, the value will be bogus.
   // In these cases, an error != 0 is returned, and we try the estimation again
   // with an increased neighborhood radius.
   int err;
   sigma =  this->EstimateSigma(idx, neighborhood,weights,pos, sigma, epsilon, err);
   
   while (err != 0)
     {
     neighborhood_radius *= 2.0;
     // Constrain the neighborhood size.  If we have reached a maximum
     // possible neighborhood size, we'll just go with that.
     if ( neighborhood_radius > this->GetMaximumNeighborhoodRadius())
       {
       sigma = this->GetMaximumNeighborhoodRadius() / this->GetNeighborhoodToSigmaRatio();
       neighborhood_radius =  this->GetMaximumNeighborhoodRadius();
       break;
       }
     else
       {
       sigma = neighborhood_radius / this->GetNeighborhoodToSigmaRatio();
       }
     
     neighborhood = system->FindNeighborhoodPoints(pos, neighborhood_radius, d);
     this->ComputeAngularWeights(pos,neighborhood,domain,weights);
     sigma = this->EstimateSigma(idx, neighborhood, weights, pos, sigma, epsilon, err);
     } // done while err
   
   // Constrain sigma to a maximum reasonable size based on the user-supplied
   // limit to neighborhood size.
   if (sigma > this->GetMaximumNeighborhoodRadius())
     {
     sigma = this->GetMaximumNeighborhoodRadius() / this->GetNeighborhoodToSigmaRatio();
     neighborhood_radius = this->GetMaximumNeighborhoodRadius();
     neighborhood = system->FindNeighborhoodPoints(pos, neighborhood_radius, d);
     this->ComputeAngularWeights(pos,neighborhood,domain,weights);
     }
 
   //  std::cout << idx <<  "\t SIGMA = " << sigma << "\t NEIGHBORHOOD SIZE = " << neighborhood.size()
   //            << "\t NEIGHBORHOOD RADIUS= " << neighborhood_radius << std::endl;
   
   // We are done with the sigma estimation step.  Cache the sigma value for
    // next time.
    m_SpatialSigmaCache->operator[](d)->operator[](idx) = sigma;
 
   //----------------------------------------------
  
    // Compute the gradients.
    double sigma2inv = 1.0 / (2.0* sigma * sigma + epsilon);
 
    VectorType r;
    VectorType gradE;
 
    for (unsigned int n = 0; n < VDimension; n++)
      {
      gradE[n] = 0.0;
      }
    
    double A = 0.0;
    for (unsigned int i = 0; i < neighborhood.size(); i++)
      {
      //    if ( neighborhood[i].Index == idx) continue;
      if (weights[i] < epsilon) continue;
      
      for (unsigned int n = 0; n < VDimension; n++)
        {
        // Note that the Neighborhood object has already filtered the
        // neighborhood for points whose normals differ by > 90 degrees.
        r[n] = pos[n] - neighborhood[i].Point[n];
        }
      
      double q = exp( -dot_product(r, r) * sigma2inv);
      A += q;
 
      for (unsigned int n = 0; n < VDimension; n++)
        {
        gradE[n] += weights[i] * r[n] * q;
        }
      }
    
    double p = 0.0;
    if (A > epsilon)
      {
      p = -1.0 / (A * sigma * sigma);
      }
 
 
    
    // TEST
    //   vnl_vector_fixed<float, VDimension> tosurf = domain->SampleGradientVnl(pos);
    //   float tosurfmag = tosurf.magnitude() + 1.0e-5;
 
    // end test
    //   float f = domain->Sample(pos);
 
    for (unsigned int n = 0; n <VDimension; n++)
      {
      gradE[n] *= p;
      // TEST
      //     gradE[n] += f * (tosurf[n] / tosurfmag);
      // end test
      }
    //   maxdt = sigma * sigma;
    maxdt = 0.5;
    
    return gradE;
 }

template<class TGradientNumericType, unsigned int VDimension>

itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::itkNewMacro ( Self )

Method for creation through the object factory.

template<class TGradientNumericType, unsigned int VDimension>

itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::itkStaticConstMacro	(	Dimension	,
		unsigned	int,
		VDimension
	)

Dimensionality of the domain of the particle system.

template<class TGradientNumericType, unsigned int VDimension>

virtual void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::ResetBuffers ( )

inlinevirtual

May be called by the solver class.

Reimplemented from itk::PSMCostFunction< VDimension >.

Definition at line 133 of file itkPSMParticleEntropyFunction.h.

   {
     m_SpatialSigmaCache->ZeroAllValues();
   }

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SetFlatCutoff ( double s )

inline

The influence of particle neighbors on each other is weighted by the angle between the surface normals at their respective positions. The "flat cutoff" parameter adjusts the degree to which surface normals must differ before the weighting kicks in. Angles below the "flat cutoff" angle, which is specified in radians, are given a weighting of 1.0. Angles above this parameter value result in a lower-weight interaction between particles. Flat cutoff is set to a default value that should work for most applications. It is not necessary to set or tune this parameter unless you would like to tweak performance on a particular dataset.

Definition at line 189 of file itkPSMParticleEntropyFunction.h.

190 { m_FlatCutoff = s; }

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SetMaximumNeighborhoodRadius ( double s )

inline

Maximum radius of the neighborhood of points that are considered in the calculation. The neighborhood is a spherical radius in 3D space. MaximumNeighborhoodRadius should be set to a value equal-to or less-than the length of the largest dimension of the image. This parameter should ALWAYS be set by an application, since this class cannot know by default the size of input images. The radius value should be specified in real-world coordinates.

Definition at line 173 of file itkPSMParticleEntropyFunction.h.

174 { m_MaximumNeighborhoodRadius = s; }

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SetMinimumNeighborhoodRadius ( double s )

inline

The minimum radius of the neighborhood of points that are considered in the entropy calculation. The neighborhood is a spherical radius in 3D space. The actual radius used in a calculation may exceed this value, but will not exceed the MaximumNeighborhoodRadius. This parameter should ALWAYS be set by an application, because it must be scaled according to the spacing in the image. A good value is typically 5x the spacing of the highest-resolution dimension (the dimension with the smallest spacing.

Definition at line 160 of file itkPSMParticleEntropyFunction.h.

161 { m_MinimumNeighborhoodRadius = s; }

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SetNeighborhoodToSigmaRatio ( double s )

inline

NeighborhoodToSigmaRatio defines the extent of a given particle's neighborhood. A particle only interacts with neighbors in this neighborhood. All other particles on the surface are ignored. The neighborhood extent is computed as a multiple of the estimated standard deviation (sigma) of the local particle positions. This parameter has a default value that should work well for most data. It is not necessary to set or tune this parameter unless you would like to tweak performance for a particular dataset.

Definition at line 203 of file itkPSMParticleEntropyFunction.h.

204 { m_NeighborhoodToSigmaRatio = s; }

template<class TGradientNumericType, unsigned int VDimension>

void itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >::SetSpatialSigmaCache ( SigmaCacheType * s )

inline

Access the cache of sigma values for each particle position. This cache is populated by registering this object as an observer of the correct particle system (see SetParticleSystem).

Definition at line 211 of file itkPSMParticleEntropyFunction.h.

212 { m_SpatialSigmaCache = s; }

The documentation for this class was generated from the following files:

C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/include/itkPSMParticleEntropyFunction.h
C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/include/itkPSMParticleEntropyFunction.hxx

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

template<class TGradientNumericType, unsigned int VDimension> class itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >

Member Typedef Documentation

Member Function Documentation

template<class TGradientNumericType, unsigned int VDimension>
class itk::PSMParticleEntropyFunction< TGradientNumericType, VDimension >