Initial commit. Taken from the Review/Statistics/ directory of ITK 3.…

…18---used to define the individual gaussian mixture model component in the Atropos segmentation.

Temporary/itkGaussianMembershipFunction.h

@@ -0,0 +1,117 @@
+/*=========================================================================
+
+  Program:   Insight Segmentation & Registration Toolkit
+  Module:    $RCSfile: itkGaussianMembershipFunction.h,v $
+  Language:  C++
+  Date:      $Date: 2009-05-02 05:43:55 $
+  Version:   $Revision: 1.1 $
+
+  Copyright (c) Insight Software Consortium. All rights reserved.
+  See ITKCopyright.txt or http://www.itk.org/HTML/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 notices for more information.
+
+=========================================================================*/
+#ifndef __itkGaussianMembershipFunction_h
+#define __itkGaussianMembershipFunction_h
+
+#include "itkArray.h"
+#include "itkMatrix.h"
+#include "itkMembershipFunctionBase.h"
+
+namespace itk
+{
+namespace Statistics
+{
+/** \class GaussianMembershipFunction
+ * \brief GaussianMembershipFunction class represents Gaussian function.
+ *
+ * This class keeps parameter to define Gaussian function and has
+ * method to return the probability density of an instance (pattern) .
+ * If the all element of the covariance matrix is zero the "usual" density
+ * calculations ignored. if the measurement vector to be evaluated is equal to
+ * the mean, then the Evaluate method will return maximum value of
+ * double and return 0 for others
+ *
+ *
+ */
+
+template <class TMeasurementVector>
+class ITK_EXPORT GaussianMembershipFunction :
+  public         MembershipFunctionBase<TMeasurementVector>
+{
+public:
+  /** Standard class typedefs */
+  typedef GaussianMembershipFunction                 Self;
+  typedef MembershipFunctionBase<TMeasurementVector> Superclass;
+  typedef SmartPointer<Self>                         Pointer;
+  typedef SmartPointer<const Self>                   ConstPointer;
+
+  /** Strandard macros */
+  itkTypeMacro(GaussianMembershipFunction, MembershipFunction);
+  itkNewMacro(Self);
+
+  /** Typedef alias for the measurement vectors */
+  typedef TMeasurementVector MeasurementVectorType;
+
+  /** Length of each measurement vector */
+  typedef typename Superclass::MeasurementVectorSizeType MeasurementVectorSizeType;
+
+  /** Type of the mean vector */
+  typedef Array<double> MeanType;
+
+  /** Type of the covariance matrix */
+  typedef VariableSizeMatrix<double> CovarianceType;
+
+  /** Set/Get the mean */
+  void SetMean( const MeanType & mean );
+
+  itkGetConstMacro( Mean, MeanType );
+
+  /** Sets the covariance matrix.
+   * Also, this function calculates inverse covariance and pre factor of
+   * Gaussian Distribution to speed up GetProbability */
+  void SetCovariance(const CovarianceType & cov);
+
+  itkGetConstMacro( Covariance, CovarianceType );
+
+  /** Gets the probability density of a measurement vector. */
+  double Evaluate(const MeasurementVectorType & measurement) const;
+
+  /** Return a copy of the current membership function */
+  Pointer Clone();
+
+protected:
+  GaussianMembershipFunction(void);
+  virtual ~GaussianMembershipFunction(void)
+  {
+  }
+
+  void PrintSelf(std::ostream& os, Indent indent) const;
+
+private:
+  MeanType       m_Mean;            // mean
+  CovarianceType m_Covariance;      // covariance matrix
+
+  // inverse covariance matrix which is automatically calculated
+  // when covariace matirx is set.  This speed up the GetProbability()
+  CovarianceType m_InverseCovariance;
+
+  // pre_factor which is automatically calculated
+  // when covariace matirx is set.  This speeds up the GetProbability()
+  double m_PreFactor;
+
+  /** if the all element of the given covarinace is zero, then this
+   * value set to true */
+  bool m_IsCovarianceZero;
+};
+} // end of namespace Statistics
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkGaussianMembershipFunction.txx"
+#endif
+
+#endif

Temporary/itkGaussianMembershipFunction.txx

@@ -0,0 +1,185 @@
+/*=========================================================================
+
+  Program:   Insight Segmentation & Registration Toolkit
+  Module:    $RCSfile: itkGaussianMembershipFunction.txx,v $
+  Language:  C++
+  Date:      $Date: 2009-08-08 15:48:18 $
+  Version:   $Revision: 1.2 $
+
+  Copyright (c) Insight Software Consortium. All rights reserved.
+  See ITKCopyright.txt or http://www.itk.org/HTML/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 notices for more information.
+
+=========================================================================*/
+#ifndef __itkGaussianMembershipFunction_txx
+#define __itkGaussianMembershipFunction_txx
+
+#include "itkGaussianMembershipFunction.h"
+
+namespace itk
+{
+namespace Statistics
+{
+template <class TMeasurementVector>
+GaussianMembershipFunction<TMeasurementVector>
+::GaussianMembershipFunction()
+{
+  m_PreFactor = 0.0;
+  m_Covariance.SetIdentity();
+}
+
+template <class TMeasurementVector>
+void
+GaussianMembershipFunction<TMeasurementVector>
+::PrintSelf(std::ostream& os, Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+
+  os << indent << "Mean: " << m_Mean << std::endl;
+  os << indent << "Covariance: " << std::endl;
+  os << m_Covariance.GetVnlMatrix();
+  os << indent << "InverseCovariance: " << std::endl;
+  os << indent << m_InverseCovariance.GetVnlMatrix();
+  os << indent << "Prefactor: " << m_PreFactor << std::endl;
+}
+
+template <class TMeasurementVector>
+void
+GaussianMembershipFunction<TMeasurementVector>
+::SetMean( const MeanType & mean )
+{
+  if( this->GetMeasurementVectorSize() )
+    {
+    MeasurementVectorTraits::Assert(mean,
+                                    this->GetMeasurementVectorSize(),
+                                    "GaussianMembershipFunction::SetMean Size of measurement vectors in \
+    the sample must the same as the size of the mean."                                                                                                          );
+    }
+  else
+    {
+    this->SetMeasurementVectorSize( mean.Size() );
+    }
+
+  if( m_Mean != mean )
+    {
+    m_Mean = mean;
+    this->Modified();
+    }
+}
+
+template <class TMeasurementVector>
+void
+GaussianMembershipFunction<TMeasurementVector>
+::SetCovariance(const CovarianceType & cov)
+{
+  // Sanity check
+  if( cov.GetVnlMatrix().rows() != cov.GetVnlMatrix().cols() )
+    {
+    itkExceptionMacro( << "Covariance matrix must be square" );
+    }
+  if( this->GetMeasurementVectorSize() )
+    {
+    if( cov.GetVnlMatrix().rows() != this->GetMeasurementVectorSize() )
+      {
+      itkExceptionMacro( << "Length of measurement vectors in the sample must be"
+                         << " the same as the size of the covariance." );
+      }
+    }
+  else
+    {
+    this->SetMeasurementVectorSize( cov.GetVnlMatrix().rows() );
+    }
+
+  m_Covariance = cov;
+
+  m_IsCovarianceZero = m_Covariance.GetVnlMatrix().is_zero();
+
+  if( !m_IsCovarianceZero )
+    {
+    // allocate the memory for m_InverseCovariance matrix
+    m_InverseCovariance.GetVnlMatrix() =
+      vnl_matrix_inverse<double>(m_Covariance.GetVnlMatrix() );
+
+    // the determinant of the covaraince matrix
+    double det = vnl_determinant(m_Covariance.GetVnlMatrix() );
+
+    // calculate coefficient C of multivariate gaussian
+    m_PreFactor = 1.0 / (vcl_sqrt(det)
+                         * vcl_pow(vcl_sqrt(2.0 * vnl_math::pi), double(this->GetMeasurementVectorSize() ) ) );
+    }
+}
+
+template <class TMeasurementVector>
+inline double
+GaussianMembershipFunction<TMeasurementVector>
+::Evaluate(const MeasurementVectorType & measurement) const
+{
+  double temp;
+
+  const MeasurementVectorSizeType measurementVectorSize =
+    this->GetMeasurementVectorSize();
+  MeanType tempVector;
+  MeasurementVectorTraits::SetLength( tempVector, measurementVectorSize );
+  MeanType tempVector2;
+  MeasurementVectorTraits::SetLength( tempVector2, measurementVectorSize );
+
+  if( !m_IsCovarianceZero )
+    {
+    // Compute |y - mean |
+    for( unsigned int i = 0; i < measurementVectorSize; i++ )
+      {
+      tempVector[i] = measurement[i] - m_Mean[i];
+      }
+    // Compute |y - mean | * inverse(cov)
+    for( unsigned int i = 0; i < measurementVectorSize; i++ )
+      {
+      temp = 0;
+      for( unsigned int j = 0; j < measurementVectorSize; j++ )
+        {
+        temp += tempVector[j] * m_InverseCovariance.GetVnlMatrix().get(j, i);
+        }
+      tempVector2[i] = temp;
+      }
+
+    // Compute |y - mean | * inverse(cov) * |y - mean|^T
+    temp = 0;
+    for( unsigned int i = 0; i < measurementVectorSize; i++ )
+      {
+      temp += tempVector2[i] * tempVector[i];
+      }
+
+    return m_PreFactor * vcl_exp(-0.5 * temp );
+    }
+  else
+    {
+    for( unsigned int i = 0; i < measurementVectorSize; i++ )
+      {
+      if( m_Mean[i] != (double) measurement[i] )
+        {
+        return 0;
+        }
+      }
+    return NumericTraits<double>::max();
+    }
+}
+
+template <class TVector>
+typename GaussianMembershipFunction<TVector>::Pointer
+GaussianMembershipFunction<TVector>
+::Clone()
+{
+  Pointer membershipFunction = GaussianMembershipFunction<TVector>::New();
+
+  membershipFunction->SetMeasurementVectorSize( this->GetMeasurementVectorSize() );
+  membershipFunction->SetMean( this->GetMean() );
+  membershipFunction->SetCovariance( this->GetCovariance() );
+
+  return membershipFunction;
+}
+} // end namespace Statistics
+} // end of namespace itk
+
+#endif