itkStochasticTractographyFilter.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 __itkStochasticTractographyFilter_h__
#define __itkStochasticTractographyFilter_h__

#include "itkImageToImageFilter.h"
#include "vnl/vnl_random.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/vnl_matrix.h"
#include "itkArray.h"
#include "itkVectorContainer.h"
#include "vnl/algo/vnl_qr.h"
#include "itkVariableLengthVector.h"
#include "StochasticTracking/itkSlowPolyLineParametricPath.h"
#include "itkSimpleFastMutexLock.h"
#include "itkRealTimeClock.h"
#include "itkDiffusionTensor3D.h"
#include <vector>

namespace itk{

typedef Image< Array< double >, 3 > ProbabilityDistributionImageType;

template< class TInputDWIImage, class TInputWhiteMatterProbabilityImage,
  class TOutputConnectivityImage >
class ITK_EXPORT StochasticTractographyFilter :
  public ImageToImageFilter< TInputDWIImage,
                    TOutputConnectivityImage >{
public:
  typedef StochasticTractographyFilter Self;
  typedef ImageToImageFilter< TInputDWIImage,
    TOutputConnectivityImage > Superclass;
  typedef SmartPointer< Self > Pointer;
  typedef SmartPointer< const Self > ConstPointer;

  itkFactorylessNewMacro(Self)
  itkCloneMacro(Self)
  itkTypeMacro( StochasticTractographyFilter,
            ImageToImageFilter )

  typedef TInputDWIImage InputDWIImageType;

  typedef TOutputConnectivityImage OutputConnectivityImageType;

  typedef TInputWhiteMatterProbabilityImage InputWhiteMatterProbabilityImageType;

  typedef SlowPolyLineParametricPath< 3 > TractType;

  typedef VectorContainer< unsigned int, typename TractType::Pointer >
    TractContainerType;

  typedef Image< DiffusionTensor3D< double >, 3 > OutputTensorImageType;

  typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > >
    GradientDirectionContainerType;

  typedef double bValueType;
  typedef VectorContainer< unsigned int, bValueType > bValueContainerType;

  typedef vnl_matrix_fixed< double, 3, 3 > MeasurementFrameType;

  typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > >
    TractOrientationContainerType;

  itkSetMacro( TotalTracts, unsigned int)
  itkGetMacro( TotalTracts, unsigned int)

  itkSetMacro( MaxTractLength, unsigned int )
  itkGetMacro( MaxTractLength, unsigned int )

  itkSetConstObjectMacro( bValues, bValueContainerType )
  itkGetConstObjectMacro( bValues, bValueContainerType )

  itkSetConstObjectMacro( Gradients, GradientDirectionContainerType )
  itkGetConstObjectMacro( Gradients, GradientDirectionContainerType )

  /* At each voxel specifies the probability of a mylinated fiber existing
     at that location.  This probability is interpreted to be the probability
     that a fiber tract passes through that region.
     */
  itkSetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType)
  itkGetInputMacro(WhiteMatterProbabilityImage, InputWhiteMatterProbabilityImageType)

  //overide the built in set input function
  //we need to create a new cache everytime we change the input image
  //but we need to preserve it when the input image is the same
  void SetPrimaryInput(DataObject *input)
  {
    typename InputDWIImageType::Pointer dwiimagePtr;
    try
    {
      dwiimagePtr = dynamic_cast<InputDWIImageType*>( input );
    }
    catch(itk::ExceptionObject &e)
    {
        MITK_INFO << "wrong image type: " << e.what();
        throw;
    }
    Superclass::SetPrimaryInput( input );
    //update the likelihood cache
    this->m_LikelihoodCachePtr = ProbabilityDistributionImageType::New();
    this->m_LikelihoodCachePtr->CopyInformation( this->GetInput() );
    this->m_LikelihoodCachePtr->SetBufferedRegion( this->GetInput()->GetBufferedRegion() );
    this->m_LikelihoodCachePtr->SetRequestedRegion( this->GetInput()->GetRequestedRegion() );
    this->m_LikelihoodCachePtr->Allocate();
    this->m_CurrentLikelihoodCacheElements = 0;
    //update the likelihoodcache mutex image
    this->m_LikelihoodCacheMutexImagePtr = LikelihoodCacheMutexImageType::New();
    this->m_LikelihoodCacheMutexImagePtr->CopyInformation( this->GetInput() );
    this->m_LikelihoodCacheMutexImagePtr->SetBufferedRegion( this->GetInput()->GetBufferedRegion() );
    this->m_LikelihoodCacheMutexImagePtr->SetRequestedRegion( this->GetInput()->GetRequestedRegion() );
    this->m_LikelihoodCacheMutexImagePtr->Allocate();
  }

  itkSetMacro( SeedIndex, typename InputDWIImageType::IndexType )
  itkGetMacro( SeedIndex, typename InputDWIImageType::IndexType )

  itkSetConstObjectMacro( SampleDirections, TractOrientationContainerType )
  itkGetConstObjectMacro( SampleDirections, TractOrientationContainerType )

  itkSetMacro( MeasurementFrame, MeasurementFrameType )
  itkGetMacro( MeasurementFrame, MeasurementFrameType )

  itkSetMacro( MaxLikelihoodCacheSize, unsigned int )
  itkGetMacro( MaxLikelihoodCacheSize, unsigned int )

  itkGetObjectMacro( OutputTractContainer, TractContainerType )

  itkGetObjectMacro( OutputTensorImage, OutputTensorImageType )

  void GenerateData();
  void GenerateTractContainerOutput( void );
  void GenerateTensorImageOutput( void );

protected:
  typedef vnl_vector_fixed< double, 7 > TensorModelParamType;

  typedef vnl_vector_fixed< double, 6 > ConstrainedModelParamType;

  typedef Image< VariableLengthVector< double >, 3 > DWIVectorImageType;

  typedef Image< Array< double >, 3 > ProbabilityDistributionImageType;

  typedef Image< SimpleFastMutexLock, 3 > LikelihoodCacheMutexImageType;

  StochasticTractographyFilter();
  virtual ~StochasticTractographyFilter();

  void LoadDefaultSampleDirections( void );

  void ProbabilisticallyInterpolate( vnl_random& randomgenerator,
    const TractType::ContinuousIndexType& cindex,
    typename InputDWIImageType::IndexType& index);

  void UpdateGradientDirections(void);
  void UpdateTensorModelFittingMatrices( void );
  void CalculateTensorModelParameters( const DWIVectorImageType::PixelType& dwivalues,
    vnl_diag_matrix<double>& W,
    TensorModelParamType& tensormodelparams);

  void CalculateConstrainedModelParameters( const TensorModelParamType& tensormodelparams,
    ConstrainedModelParamType& constrainedmodelparams);

  void CalculateNoiseFreeDWIFromConstrainedModel( const ConstrainedModelParamType& constrainedmodelparams,
    DWIVectorImageType::PixelType& noisefreedwi);

  void CalculateResidualVariance( const DWIVectorImageType::PixelType& noisydwi,
    const DWIVectorImageType::PixelType& noisefreedwi,
    const vnl_diag_matrix< double >& W,
    const unsigned int numberofparameters,
    double& residualvariance);

  void CalculateLikelihood( const DWIVectorImageType::PixelType &dwipixel,
    TractOrientationContainerType::ConstPointer orientations,
    ProbabilityDistributionImageType::PixelType& likelihood);

  void CalculatePrior( TractOrientationContainerType::Element v_prev,
    TractOrientationContainerType::ConstPointer orientations,
    ProbabilityDistributionImageType::PixelType& prior );

  void CalculatePosterior( const ProbabilityDistributionImageType::PixelType& likelihood,
    const ProbabilityDistributionImageType::PixelType& prior,
    ProbabilityDistributionImageType::PixelType& posterior);

  void SampleTractOrientation( vnl_random& randomgenerator,
    const ProbabilityDistributionImageType::PixelType& posterior,
    TractOrientationContainerType::ConstPointer orientations,
    TractOrientationContainerType::Element& choosendirection );

  void StochasticTractGeneration( typename InputDWIImageType::ConstPointer dwiimagePtr,
    typename InputWhiteMatterProbabilityImageType::ConstPointer maskimagePtr,
    typename InputDWIImageType::IndexType seedindex,
    unsigned long randomseed,
    TractType::Pointer tract );

  static ITK_THREAD_RETURN_TYPE StochasticTractGenerationCallback( void *arg );

  struct StochasticTractGenerationCallbackStruct{
    Pointer Filter;
  };

  ProbabilityDistributionImageType::PixelType&
    AccessLikelihoodCache( typename InputDWIImageType::IndexType index );
  bool DelegateTract(unsigned long& randomseed);
  void TractContainerToConnectivityMap(TractContainerType::Pointer tractcontainer);
  void StoreTract(TractType::Pointer tract);
  bool FiberExistenceTest( vnl_random& randomgenerator,
    typename InputWhiteMatterProbabilityImageType::ConstPointer wmpimage,
    typename InputWhiteMatterProbabilityImageType::IndexType index );



  MeasurementFrameType m_MeasurementFrame;

  LikelihoodCacheMutexImageType::Pointer m_LikelihoodCacheMutexImagePtr;
  unsigned int m_TotalTracts;
  unsigned int m_MaxTractLength;
  GradientDirectionContainerType::ConstPointer m_Gradients;

  GradientDirectionContainerType::Pointer m_TransformedGradients;
  bValueContainerType::ConstPointer m_bValues;

  typename InputDWIImageType::IndexType m_SeedIndex;
  TractOrientationContainerType::ConstPointer m_SampleDirections;

    //these will be the same for every pixel in the image so
  //go ahead and do a QR decomposition to optimize the
  //LS fitting process for estimating the weighing matrix W
  //in this case we solve instead:
  //R*Beta = Q'logPhi
  vnl_matrix< double >* m_A;
  vnl_qr< double >* m_Aqr;
  ProbabilityDistributionImageType::Pointer m_LikelihoodCachePtr;

  unsigned long m_MaxLikelihoodCacheSize;   //in Megabytes
  unsigned long m_MaxLikelihoodCacheElements;  //in Elements (Voxels)
  unsigned long m_CurrentLikelihoodCacheElements;
  SimpleFastMutexLock m_LikelihoodCacheMutex;

  RealTimeClock::Pointer m_ClockPtr;
  unsigned int m_TotalDelegatedTracts;
  SimpleFastMutexLock m_TotalDelegatedTractsMutex;

  //unsigned long m_RandomSeed;
  SimpleFastMutexLock m_OutputImageMutex;
  TractContainerType::Pointer m_OutputTractContainer;
  SimpleFastMutexLock m_OutputTractContainerMutex;

  OutputTensorImageType::Pointer m_OutputTensorImage;
  vnl_random m_RandomGenerator;
};

}

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkStochasticTractographyFilter.txx"
#include "StochasticTracking/itkStochasticTractographyFilter_SD.txx"
#endif

#endif