itkAnalyticalDiffusionQballReconstructionImageFilter.h

Go to the documentation of this file.

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

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.

This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.

See LICENSE.txt or http://www.mitk.org for details.

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

#ifndef __itkAnalyticalDiffusionQballReconstructionImageFilter_h_
#define __itkAnalyticalDiffusionQballReconstructionImageFilter_h_

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


namespace itk{
template< class TReferenceImagePixelType,
          class TGradientImagePixelType,
          class TOdfPixelType,
          int NOrderL,
          int NrOdfDirections>
class AnalyticalDiffusionQballReconstructionImageFilter :
        public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >,
        Image< Vector< TOdfPixelType, 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, NrOdfDirections >, 3 > > Superclass;

    itkFactorylessNewMacro(Self)
    itkCloneMacro(Self)

    itkTypeMacro(AnalyticalDiffusionQballReconstructionImageFilter, ImageToImageFilter)

    typedef TReferenceImagePixelType                 ReferencePixelType;

    typedef TGradientImagePixelType                  GradientPixelType;

    typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType;

    typedef TOdfPixelType                            BZeroPixelType;

    typedef typename Superclass::InputImageType       ReferenceImageType;

    typedef Image< OdfPixelType, 3 >                  OdfImageType;

    typedef OdfImageType                              OutputImageType;

    typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType;

    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( const 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 void Cart2Sph(double x, double y, double z, double* cart);
    static double Yj(int m, int k, double theta, double phi, bool useMRtrixBasis = false);
    double Legendre0(int l);

    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 )

    typedef Image<float, 3> FloatImageType;
    itkGetMacro( BZeroImage, typename BZeroImageType::Pointer)
    itkGetMacro( ODFSumImage, typename FloatImageType::Pointer)
    itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer)

    itkSetMacro( BValue, TOdfPixelType)
#ifdef GetBValue
#undef GetBValue
#endif
    itkGetConstReferenceMacro( BValue, TOdfPixelType)
    itkSetMacro( Lambda, double )
    itkGetMacro( Lambda, double )

    itkSetMacro( UseMrtrixBasis, bool )

#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, ThreadIdType);

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;
    vnl_matrix<double>*                               m_B_t;
    vnl_vector<double>*                               m_LP;
    FloatImageType::Pointer                           m_ODFSumImage;
    typename CoefficientImageType::Pointer            m_CoefficientImage;
    TOdfPixelType                                     m_Delta1;
    TOdfPixelType                                     m_Delta2;
    bool                                              m_UseMrtrixBasis;
};

}

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

#endif //__itkAnalyticalDiffusionQballReconstructionImageFilter_h_