itkAnalyticalDiffusionQballReconstructionImageFilter.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2006-03-27 17:01:06 $
  Version:   $Revision: 1.12 $

  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 __itkAnalyticalDiffusionQballReconstructionImageFilter_h_
#define __itkAnalyticalDiffusionQballReconstructionImageFilter_h_

#include "itkImageToImageFilter.h"
//#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/vnl_matrix.h"
#include "vnl/algo/vnl_svd.h"
#include "itkVectorContainer.h"
#include "itkVectorImage.h"
#include "QuadProg.h"

namespace itk{
template< class TReferenceImagePixelType, 
          class TGradientImagePixelType,
          class TOdfPixelType,
          int NOrderL,
          int NrOdfDirections>
class AnalyticalDiffusionQballReconstructionImageFilter :
  public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, 
                             Image< Vector< TOdfPixelType, /*(NOrderL*NOrderL + NOrderL + 2.0)/2.0 + NOrderL*/NrOdfDirections >, 3 > >
{

public:

  enum Normalization {
    QBAR_STANDARD,
    QBAR_B_ZERO_B_VALUE,
    QBAR_B_ZERO,
    QBAR_NONE,
    QBAR_ADC_ONLY,
    QBAR_RAW_SIGNAL,
    QBAR_SOLID_ANGLE,
    QBAR_NONNEG_SOLID_ANGLE
  };

  typedef AnalyticalDiffusionQballReconstructionImageFilter Self;
  typedef SmartPointer<Self>                      Pointer;
  typedef SmartPointer<const Self>                ConstPointer;
  typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, 
          Image< Vector< TOdfPixelType, /*NOrderL*NOrderL + NOrderL + 2.0)/2.0 + NOrderL*/NrOdfDirections >, 3 > >
                          Superclass;
  
  itkNewMacro(Self);  

  itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, 
                                                   ImageToImageFilter);
 
  typedef TReferenceImagePixelType                 ReferencePixelType;

  typedef TGradientImagePixelType                  GradientPixelType;

  typedef Vector< TOdfPixelType, /*NOrderL*NOrderL + NOrderL + 2.0)/2.0 + NOrderL*/NrOdfDirections >
                                                   OdfPixelType;

  typedef typename Superclass::InputImageType      ReferenceImageType;
  
  typedef Image< OdfPixelType, 3 >                 OdfImageType;

  typedef OdfImageType                             OutputImageType;

  typedef TOdfPixelType                            BZeroPixelType;

  typedef Image< BZeroPixelType, 3 >      BZeroImageType;

  typedef typename Superclass::OutputImageRegionType
                                                   OutputImageRegionType;

  typedef Image< GradientPixelType, 3 >            GradientImageType;

  typedef VectorImage< GradientPixelType, 3 >      GradientImagesType;

  typedef vnl_matrix< TOdfPixelType >*
                                                   OdfReconstructionMatrixType;
  
  typedef vnl_matrix< double >                     CoefficientMatrixType;

  typedef vnl_vector_fixed< double, 3 >            GradientDirectionType;

  typedef VectorContainer< unsigned int, 
          GradientDirectionType >                  GradientDirectionContainerType;
  
  void SetGradientImage( GradientDirectionContainerType *, 
                                             const GradientImagesType *image);
  
  virtual ReferenceImageType * GetReferenceImage() 
  { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); }

  virtual GradientDirectionType GetGradientDirection( unsigned int idx) const
    {
    if( idx >= m_NumberOfGradientDirections )
      {
      itkExceptionMacro( << "Gradient direction " << idx << "does not exist" );
      }
    return m_GradientDirectionContainer->ElementAt( idx+1 );
    }

  static void tofile2(vnl_matrix<double> *A, std::string fname);
  static double factorial(int number);
  static void Cart2Sph(double x, double y, double z, double* cart);
  static double legendre0(int l);
  static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart);
  static double Yj(int m, int k, double theta, double phi);

  OdfPixelType Normalize(OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 );

  vnl_vector<TOdfPixelType> PreNormalize( vnl_vector<TOdfPixelType> vec, typename NumericTraits<ReferencePixelType>::AccumulateType b0  );

  itkSetMacro( Threshold, ReferencePixelType );
  itkGetMacro( Threshold, ReferencePixelType );

  itkSetMacro( NormalizationMethod, Normalization);
  itkGetMacro( NormalizationMethod, Normalization );

  itkGetMacro( BZeroImage, typename BZeroImageType::Pointer);

  itkSetMacro( BValue, TOdfPixelType);
#ifdef GetBValue
#undef GetBValue
#endif
  itkGetConstReferenceMacro( BValue, TOdfPixelType);

  itkSetMacro( Lambda, double );
  itkGetMacro( Lambda, double );

#ifdef ITK_USE_CONCEPT_CHECKING

  itkConceptMacro(ReferenceEqualityComparableCheck,
    (Concept::EqualityComparable<ReferencePixelType>));
  itkConceptMacro(TensorEqualityComparableCheck,
    (Concept::EqualityComparable<OdfPixelType>));
  itkConceptMacro(GradientConvertibleToDoubleCheck,
    (Concept::Convertible<GradientPixelType, double>));
  itkConceptMacro(DoubleConvertibleToTensorCheck,
    (Concept::Convertible<double, OdfPixelType>));
  itkConceptMacro(GradientReferenceAdditiveOperatorsCheck,
    (Concept::AdditiveOperators<GradientPixelType, GradientPixelType,
                                ReferencePixelType>));
  itkConceptMacro(ReferenceOStreamWritableCheck,
    (Concept::OStreamWritable<ReferencePixelType>));
  itkConceptMacro(TensorOStreamWritableCheck,
    (Concept::OStreamWritable<OdfPixelType>));
#endif

protected:
  AnalyticalDiffusionQballReconstructionImageFilter();
  ~AnalyticalDiffusionQballReconstructionImageFilter() {};
  void PrintSelf(std::ostream& os, Indent indent) const;

  void ComputeReconstructionMatrix();
  
  void BeforeThreadedGenerateData();
  void ThreadedGenerateData( const 
      OutputImageRegionType &outputRegionForThread, int);
  
private:
  
  OdfReconstructionMatrixType                       m_ReconstructionMatrix;
  
  OdfReconstructionMatrixType                       m_CoeffReconstructionMatrix;

  OdfReconstructionMatrixType                       m_SphericalHarmonicBasisMatrix;

  GradientDirectionContainerType::Pointer           m_GradientDirectionContainer;

  unsigned int                                      m_NumberOfGradientDirections;

  unsigned int                                      m_NumberOfBaselineImages;

  ReferencePixelType                                m_Threshold;

  TOdfPixelType                                     m_BValue;

  typename BZeroImageType::Pointer                  m_BZeroImage;

  double                                            m_Lambda;

  bool                                              m_DirectionsDuplicated;

  Normalization                                     m_NormalizationMethod;

  int                                               m_NumberCoefficients;

  QuadProgPP::Matrix<double>                        m_G, m_CE, m_CI;
  QuadProgPP::Vector<double>                        m_g0, m_ce0, m_ci0, m_x;
  vnl_matrix<double>*                               m_B_t;
  vnl_vector<double>*                               m_LP;
  
};

}

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAnalyticalDiffusionQballReconstructionImageFilter.cpp"
#endif

#endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_