itkPSMEntropyModelFilter.h

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef __itkPSMEntropyModelFilter_h
#define __itkPSMEntropyModelFilter_h

#include "itkInPlaceImageFilter.h"
#include "itkPSMParticleSystem.h"
#include "itkPSMImplicitSurfaceDomain.h"
#include "itkPSMShapeEntropyFunction.h"
#include "itkPSMParticleEntropyFunction.h"
#include "itkPSMGradientDescentOptimizer.h"
#include "itkPSMContainerArrayAttribute.h"
//#include "itkPSMMeanCurvatureAttribute.h"
#include "itkPSMSurfaceNeighborhood.h"
#include "itkPSMTwoCostFunction.h"
#include "itkCommand.h"

namespace itk
{
  
template <class TImage, class TShapeMatrix = PSMShapeMatrixAttribute<double, TImage::ImageDimension> >
class ITK_EXPORT PSMEntropyModelFilter
  : public InPlaceImageFilter<TImage,TImage> 
{
 public:
  typedef PSMEntropyModelFilter  Self;
  typedef InPlaceImageFilter<TImage,TImage>  Superclass;
  typedef SmartPointer<Self>   Pointer;
  typedef SmartPointer<const Self>  ConstPointer;
  
  itkStaticConstMacro(Dimension, unsigned int, TImage::ImageDimension);
  
  typedef PSMParticleSystem<Dimension> ParticleSystemType;
  
  typedef TShapeMatrix ShapeMatrixType;
  typedef typename ShapeMatrixType::Pointer ShapeMatrixPointerType;
  
  itkNewMacro(Self);
  
  itkTypeMacro(PSMEntropyModelFilter, InPlaceImageFilter);
  
  typedef TImage ImageType;
  
  typedef typename ImageType::PointType PointType;
  
  typedef PSMGradientDescentOptimizer<typename ImageType::PixelType, Dimension> OptimizerType;

  struct CuttingPlaneType
  {
    vnl_vector_fixed<double,Dimension> x;
    vnl_vector_fixed<double,Dimension> y;
    vnl_vector_fixed<double,Dimension> z;
    bool valid;

    CuttingPlaneType()
    { valid = false; }
    
  };
  
  void SetInput(const std::string& s, itk::DataObject *o)
    {
      Superclass::SetInput(s,o);
    }
  void SetInput(const TImage *image)
  { 
    this->SetInput(0, image);  
  }  
  void SetInput( unsigned int index, const TImage * image ) 
    {
      if (this->GetNumberOfInputs() < index+1)
        {
          this->SetNumberOfRequiredInputs(index+1);
        }
      
      this->ProcessObject::SetNthInput(index, const_cast< TImage *>( image ) );
    }
  
  void SetNumberOfScales(unsigned int n)
    {
      m_NumberOfScales = n;
      if (m_Tolerances.size() < n)
        {
          m_Tolerances.resize(n);
        }
      if (m_MaximumNumberOfIterations.size() < n)
        {
          m_MaximumNumberOfIterations.resize(n);
        }
      if (m_RegularizationInitial.size() < n)
        {
          m_RegularizationInitial.resize(n);
        }
      if (m_RegularizationFinal.size() < n)
        {
          m_RegularizationFinal.resize(n);
        }
      if (m_RegularizationDecaySpan.size() < n)
        {
          m_RegularizationDecaySpan.resize(n);
        }

    }
  itkGetMacro(NumberOfScales, unsigned int);
  itkSetMacro(CurrentScale, unsigned int);
  itkGetMacro(CurrentScale, unsigned int);
  
  //  void SetInputModel(const ParticleModel &model);
  
  void SetInputCorrespondencePoints(unsigned int index, const std::vector<PointType> &corr)
    {
      if (m_InputCorrespondencePoints.size() < (index+1))
        {
          m_InputCorrespondencePoints.resize(index+1);
        }
      m_InputCorrespondencePoints[index] = corr;
    }
  
  itkGetObjectMacro(ParticleSystem, ParticleSystemType);
  itkGetConstObjectMacro(ParticleSystem, ParticleSystemType);
  
  PSMParticleEntropyFunction<typename ImageType::PixelType, Dimension>
    *GetParticleEntropyFunction()
  {
    return m_ParticleEntropyFunction;
  }
  
  PSMShapeEntropyFunction<Dimension> *ShapeEntropyFunction()
  {
    return m_ShapeEntropyFunction;
  }
  
  itkGetObjectMacro(Optimizer, OptimizerType);
  itkGetConstObjectMacro(Optimizer, OptimizerType);
  
  void SetMaximumNumberOfIterations(const std::vector<unsigned int> &n)
    {
      m_MaximumNumberOfIterations = n;
    }
  void SetMaximumNumberOfIterations(unsigned int n, unsigned int i = 0)
    {
      // Resize the tolerances list if necessary
      if (m_MaximumNumberOfIterations.size() < i+1)
        {
          itkExceptionMacro("Specified scale does not exist.  Call SetNumberOfScales before SetMaximumNumberOfIterations");
        }
      m_MaximumNumberOfIterations[i] = n;
      
      if (m_CurrentScale == i)
        {
          m_Optimizer->SetMaximumNumberOfIterations(n);
        }
      
    }
  unsigned int GetMaximumNumberOfIterations(unsigned int i = 0) const
  {
    if (m_MaximumNumberOfIterations.size() < i+1)
      {
        itkExceptionMacro("Specified scale does not exist.");
      }
    return m_MaximumNumberOfIterations[i];
  }
  
  unsigned int GetNumberOfElapsedIterations() const
  {
    return m_Optimizer->GetNumberOfIterations();
  }

  itkGetMacro(Initialized, bool);

  itkGetMacro(Initializing, bool);
  itkSetMacro(Initializing, bool);

  void SetShapeEntropyWeighting(double w) 
  {
    m_CostFunction->SetRelativeGradientScaling(w);
  }
  double GetShapeEntropyWeighting() const
  {
    return m_CostFunction->GetRelativeGradientScaling();
  }

 /* /\** This method simply calls the appropriate method in */
 /*     PSMParticleEntropyFunction.  It is exposed here for */
 /*     convenience. 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. *\/ */
 /*  void SetMinimumNeighborhoodRadius( double s) */
 /*  { this->m_ParticleEntropyFunction->SetMinimumNeighborhoodRadius(s); } */
 /*  double GetMinimumNeighborhoodRadius() const */
 /*  { return this->m_ParticleEntropyFunction->GetMinimumNeighborhoodRadius(); } */

 /*  /\** 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. *\/ */
 /*  void SetMaximumNeighborhoodRadius( double s) */
 /*  {  this->m_ParticleEntropyFunction->SetMaximumNeighborhoodRadius(s); } */
 /*  double GetMaximumNeighborhoodRadius() const */
 /*  { return  this->m_ParticleEntropyFunction->GetMaximumNeighborhoodRadius(); } */

   const std::vector<double> &GetShapePCAVariances() const
   { 
     return m_ShapeEntropyFunction->GetShapePCAVariances(); 
   }
  
   void SetRegularizationInitial(const std::vector<double> &v)
     {
       if (v.size() != m_RegularizationInitial.size()) 
         {
           itkExceptionMacro("Number of variables does not match number of scales.  Call SetNumberOfScales first.");
         }
       m_RegularizationInitial = v;
     }
   void SetRegularizationInitial(double t, unsigned int i=0)
     {
       // Resize the tolerances list if necessary
       if (m_RegularizationInitial.size() < i+1)
         {
           itkExceptionMacro("Specified scale does not exist. Call SetNumberOfScales before SetTolerance");
         }
       m_RegularizationInitial[i] = t;
       
       if (m_CurrentScale == i)
         {
           m_Optimizer->SetTolerance(t);
         }       
     }
   double GetRegularizationInitial(unsigned int i=0) const
     {
       if (m_RegularizationInitial.size() < i+1)
         {
           itkExceptionMacro("Specified scale does not exist");
         }
       else return m_RegularizationInitial[i];
     }
   void SetRegularizationFinal(const std::vector<double> &v)
     {
       if (v.size() != m_RegularizationFinal.size()) 
         {
           itkExceptionMacro("Number of variables does not match number of scales.  Call SetNumberOfScales first.");
         }
       m_RegularizationFinal = v;
     }
   void SetRegularizationFinal(double t, unsigned int i=0)
     {
       // Resize the tolerances list if necessary
       if (m_RegularizationFinal.size() < i+1)
         {
           itkExceptionMacro("Specified scale does not exist. Call SetNumberOfScales before SetTolerance");
         }
       m_RegularizationFinal[i] = t;
       
       if (m_CurrentScale == i)
         {
           m_Optimizer->SetTolerance(t);
         }       
     }
   double GetRegularizationFinal(unsigned int i=0) const
   {
     if (m_RegularizationFinal.size() < i+1)
       {
         itkExceptionMacro("Specified scale does not exist");
       }
     else return m_RegularizationFinal[i];
   }
   void SetRegularizationDecaySpan(const std::vector<double> &v)
     {
       if (v.size() != m_RegularizationDecaySpan.size()) 
         {
           itkExceptionMacro("Number of variables does not match number of scales.  Call SetNumberOfScales first.");
         }
       m_RegularizationDecaySpan = v;
     }
   void SetRegularizationDecaySpan(double t, unsigned int i=0)
     {
       // Resize the tolerances list if necessary
       if (m_RegularizationDecaySpan.size() < i+1)
         {
           itkExceptionMacro("Specified scale does not exist. Call SetNumberOfScales before SetTolerance");
         }
       m_RegularizationDecaySpan[i] = t;
       
       if (m_CurrentScale == i)
         {
           m_Optimizer->SetTolerance(t);
         }       
     }
   double GetRegularizationDecaySpan(unsigned int i=0) const
   {
     if (m_RegularizationDecaySpan.size() < i+1)
       {
         itkExceptionMacro("Specified scale does not exist");
       }
     else return m_RegularizationDecaySpan[i];
   }

   void SetRegularizationConstant(double r)
   {
     m_ShapeEntropyFunction->SetMinimumVariance(r);
   }
   void SetMinimumVariance(double r)
   {
     m_ShapeEntropyFunction->SetMinimumVariance(r);
   }
   double GetRegularizationConstant() const
   {
     return m_ShapeEntropyFunction->GetMinimumVariance();
   }
   double GetMinimumVariance() const
   {
     return m_ShapeEntropyFunction->GetMinimumVariance();
   }
   
   void SetTolerance(const std::vector<double> &v)
     {
       if (v.size() != m_Tolerances.size()) 
         {
           itkExceptionMacro("Number of variables does not match number of scales.  Call SetNumberOfScales first.");
         }
       m_Tolerances = v;
     }
   void SetTolerance(double t, unsigned int i = 0)
     {
       // Resize the tolerances list if necessary
       if (m_Tolerances.size() < i+1)
         {
           itkExceptionMacro("Specified scale does not exist. Call SetNumberOfScales before SetTolerance");
         }
       m_Tolerances[i] = t;
       
       if (m_CurrentScale == i)
         {
           m_Optimizer->SetTolerance(t);
         }
     }
   double GetTolerance(unsigned int i = 0) const
     {
       if (m_Tolerances.size() > i)
         {
           return m_Tolerances[i];
         }
       else 
         {
           itkExceptionMacro("Specified scale does not exist.");
         }
     }
   
   itkGetMacro(ShapeMatrix, ShapeMatrixPointerType);
   const ShapeMatrixType *GetShapeMatrix() const
   {
     return m_ShapeMatrix;
   }
   
   bool CreateSingleCorrespondence();

   void SetDomainName(const std::string &s, unsigned int i);

   unsigned int GetDomainIndexByName(const std::string &name);

   void AddCuttingPlane(const vnl_vector_fixed<double,3> &x,
               const vnl_vector_fixed<double,3> &y,
               const vnl_vector_fixed<double,3> &z,
               unsigned int domain);
   void AddCuttingPlane(const vnl_vector_fixed<double,3> &x,
               const vnl_vector_fixed<double,3> &y,
               const vnl_vector_fixed<double,3> &z,
               const std::string &s)
   {
     this->AddCuttingPlane(x,y,z,this->GetDomainIndexByName(s));
   }

 protected:
   PSMEntropyModelFilter();
   virtual ~PSMEntropyModelFilter() {};
   
   void PrintSelf(std::ostream& os, Indent indent) const
     {
       Superclass::PrintSelf(os, indent);
     }
   
   void GenerateData();
   
   virtual void AllocateDataCaches();
   
   virtual void AllocateWorkingImages();
   
   virtual void AllocateDomainsAndNeighborhoods();
   
   void InitializeCorrespondences();
   
   itkSetMacro(Initialized, bool);
   
   void OptimizerIterateCallback(Object *, const EventObject &);

   typename PSMImplicitSurfaceDomain<
     typename ImageType::PixelType, Dimension>::Pointer
     &GetDomainByName(const std::string &name)
   {
     return this->GetDomain(this->GetDomainIndexByName(name));
   }

   typename PSMImplicitSurfaceDomain<
     typename ImageType::PixelType, Dimension>::Pointer
     &GetDomain(unsigned int i)
   {
     if (i >= m_DomainList.size())
       { itkExceptionMacro("Requested domain not available"); }
     return m_DomainList[i];
   }
   
 private:
   PSMEntropyModelFilter(const Self&); //purposely not implemented
   void operator=(const Self&);        //purposely not implemented
   
   std::vector<typename TImage::Pointer> m_WorkingImages;
   typename OptimizerType::Pointer       m_Optimizer;
   
   std::vector<unsigned int> m_MaximumNumberOfIterations;
   
   std::vector<double> m_RegularizationInitial;
   std::vector<double> m_RegularizationFinal;
   std::vector<double> m_RegularizationDecaySpan;
   
   bool m_Initializing;
   
  bool m_Initialized;
  
  virtual void InitializeOptimizationFunctions();
  
  typename PSMContainerArrayAttribute<double, Dimension>::Pointer m_SigmaCache;
  //  typename PSMMeanCurvatureAttribute<typename ImageType::PixelType, Dimension>
  //    ::Pointer m_MeanCurvatureCache;
  
  typename PSMShapeEntropyFunction<Dimension>::Pointer m_ShapeEntropyFunction;
  
  typename PSMParticleEntropyFunction<typename TImage::PixelType, Dimension>::Pointer
    m_ParticleEntropyFunction;
  
  typename PSMTwoCostFunction<Dimension>::Pointer m_CostFunction;
  
  ShapeMatrixPointerType m_ShapeMatrix;
  
  typename itk::MemberCommand< PSMEntropyModelFilter<TImage,TShapeMatrix> >
    ::Pointer m_IterateCallback;
  
  typename ParticleSystemType::Pointer m_ParticleSystem;
  
  std::vector<std::vector<PointType> > m_InputCorrespondencePoints;

  std::vector<std::string> m_DomainListNames;
 
  std::vector<typename PSMImplicitSurfaceDomain<
  typename ImageType::PixelType, Dimension>::Pointer> m_DomainList;
 
  std::vector<typename PSMSurfaceNeighborhood<ImageType>::Pointer> m_NeighborhoodList;
  
  unsigned int m_NumberOfScales;
  
  unsigned int m_CurrentScale;
  
  std::vector<double> m_Tolerances;

  std::vector<CuttingPlaneType> m_CuttingPlanes;
};
  
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkPSMEntropyModelFilter.hxx"
#endif

#endif