diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h index e9340d04e3..95ff207baa 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h @@ -1,121 +1,122 @@ /*=================================================================== 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 _MITK_AstroStickModel_H #define _MITK_AstroStickModel_H #include #include namespace mitk { /** * \brief Generates the diffusion signal using a collection of idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class AstroStickModel : public DiffusionSignalModel< ScalarType > { public: - AstroStickModel(); - template< class OtherType >AstroStickModel(AstroStickModel* model) + AstroStickModel(); + template< class OtherType >AstroStickModel(AstroStickModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); + + this->m_BValue = model->GetBvalue(); + this->m_Diffusivity = model->GetDiffusivity(); + this->m_Sticks = model->GetSticks(); + this->m_NumSticks = model->GetNumSticks(); + this->m_RandomizeSticks = model->GetRandomizeSticks(); + } + ~AstroStickModel(); + + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + + void SetRandomizeSticks(bool randomize=true){ m_RandomizeSticks=randomize; } ///< Random stick configuration in each voxel + bool GetRandomizeSticks() { return m_RandomizeSticks; } + + void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant + double GetDiffusivity() { return m_Diffusivity; } + + void SetNumSticks(unsigned int order) + { + vnl_matrix sticks; + switch (order) { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); - - this->m_BValue = model->GetBvalue(); - this->m_Diffusivity = model->GetDiffusivity(); - this->m_Sticks = model->GetSticks(); - this->m_NumSticks = model->GetNumSticks(); - this->m_RandomizeSticks = model->GetRandomizeSticks(); + case 1: + m_NumSticks = 12; + sticks = itk::PointShell<12, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; + case 2: + m_NumSticks = 42; + sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; + case 3: + m_NumSticks = 92; + sticks = itk::PointShell<92, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; + case 4: + m_NumSticks = 162; + sticks = itk::PointShell<162, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; + case 5: + m_NumSticks = 252; + sticks = itk::PointShell<252, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; + default: + m_NumSticks = 42; + sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); + break; } - ~AstroStickModel(); - - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); - - void SetRandomizeSticks(bool randomize=true){ m_RandomizeSticks=randomize; } ///< Random stick configuration in each voxel - bool GetRandomizeSticks() { return m_RandomizeSticks; } - - void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant - double GetDiffusivity() { return m_Diffusivity; } - - void SetNumSticks(unsigned int order) + for (int i=0; i sticks; - switch (order) - { - case 1: - m_NumSticks = 12; - sticks = itk::PointShell<12, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - case 2: - m_NumSticks = 42; - sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - case 3: - m_NumSticks = 92; - sticks = itk::PointShell<92, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - case 4: - m_NumSticks = 162; - sticks = itk::PointShell<162, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - case 5: - m_NumSticks = 252; - sticks = itk::PointShell<252, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - default: - m_NumSticks = 42; - sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); - break; - } - for (int i=0; i namespace mitk { /** * \brief Generates direction independent diffusion measurement employing a scalar diffusion constant d: e^(-bd) * */ template< class ScalarType = double > class BallModel : public DiffusionSignalModel< ScalarType > { public: - BallModel(); - template< class OtherType >BallModel(BallModel* model) - { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); + BallModel(); + template< class OtherType >BallModel(BallModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); - this->m_BValue = model->GetBvalue(); - this->m_Diffusivity = model->GetDiffusivity(); - } - ~BallModel(); + this->m_BValue = model->GetBvalue(); + this->m_Diffusivity = model->GetDiffusivity(); + } + ~BallModel(); - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); - void SetDiffusivity(double D) { m_Diffusivity = D; } - double GetDiffusivity() { return m_Diffusivity; } + void SetDiffusivity(double D) { m_Diffusivity = D; } + double GetDiffusivity() { return m_Diffusivity; } protected: - double m_Diffusivity; ///< Scalar diffusion constant + double m_Diffusivity; ///< Scalar diffusion constant }; } #include "mitkBallModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h index 721ef14ad5..004fb2f417 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h @@ -1,62 +1,63 @@ /*=================================================================== 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 _MITK_DotModel_H #define _MITK_DotModel_H #include namespace mitk { /** * \brief Generates constant direction independent signal. * */ template< class ScalarType = double > class DotModel : public DiffusionSignalModel< ScalarType > { public: - DotModel(); - template< class OtherType >DotModel(DotModel* model) - { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); - } - ~DotModel(); - - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + DotModel(); + template< class OtherType >DotModel(DotModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); + } + ~DotModel(); + + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); protected: }; } #include "mitkDotModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h index 8d7f80c9fe..8f321c3867 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h @@ -1,109 +1,110 @@ /*=================================================================== 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 _MITK_RawShModel_H #define _MITK_RawShModel_H #include #include #include namespace mitk { /** * \brief The spherical harmonic representation of a prototype diffusion weighted MR signal is used to obtain the direction dependent signal. * */ template< class ScalarType = double > class RawShModel : public DiffusionSignalModel< ScalarType > { public: - RawShModel(); - template< class OtherType >RawShModel(RawShModel* model) - { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); - - this->m_AdcRange = model->GetAdcRange(); - this->m_FaRange = model->GetFaRange(); - this->m_ShCoefficients = model->GetShCoefficients(); - this->m_B0Signal = model->GetB0Signals(); - this->m_ShOrder = model->GetShOrder(); - this->m_ModelIndex = model->GetModelIndex(); - this->m_MaxNumKernels = model->GetMaxNumKernels(); - } - ~RawShModel(); - - typedef itk::Image< double, 3 > ItkDoubleImageType; - typedef itk::Image< unsigned char, 3 > ItkUcharImageType; - typedef itk::Image< itk::DiffusionTensor3D< double >, 3 > TensorImageType; - typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - typedef itk::Matrix< double, 3, 3 > MatrixType; - - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); - - bool SetShCoefficients(vnl_vector< double > shCoefficients, double b0); - vnl_matrix SetFiberDirection(GradientType& fiberDirection); - void SetFaRange(double min, double max){ m_FaRange.first = min; m_FaRange.second = max; } - void SetAdcRange(double min, double max){ m_AdcRange.first = min; m_AdcRange.second = max; } - void SetMaxNumKernels(unsigned int max){ m_MaxNumKernels = max; } - unsigned int GetNumberOfKernels(); - std::pair< double, double > GetFaRange(){ return m_FaRange; } - std::pair< double, double > GetAdcRange(){ return m_AdcRange; } - unsigned int GetMaxNumKernels(){ return m_MaxNumKernels; } - void Clear(); - - std::vector< vnl_vector< double > > GetShCoefficients(){ return m_ShCoefficients; } - std::vector< double > GetB0Signals(){ return m_B0Signal; } - unsigned int GetShOrder(){ return m_ShOrder; } - int GetModelIndex(){ return m_ModelIndex; } - - double GetBaselineSignal(int index){ return m_B0Signal.at(index); } - vnl_vector< double > GetCoefficients(int listIndex){ return m_ShCoefficients.at(listIndex); } - - bool SampleKernels(Image::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage=nullptr, QballFilterType::CoefficientImageType::Pointer coeffImage=nullptr, ItkDoubleImageType::Pointer adcImage=nullptr); + RawShModel(); + template< class OtherType >RawShModel(RawShModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); + + this->m_AdcRange = model->GetAdcRange(); + this->m_FaRange = model->GetFaRange(); + this->m_ShCoefficients = model->GetShCoefficients(); + this->m_B0Signal = model->GetB0Signals(); + this->m_ShOrder = model->GetShOrder(); + this->m_ModelIndex = model->GetModelIndex(); + this->m_MaxNumKernels = model->GetMaxNumKernels(); + } + ~RawShModel(); + + typedef itk::Image< double, 3 > ItkDoubleImageType; + typedef itk::Image< unsigned char, 3 > ItkUcharImageType; + typedef itk::Image< itk::DiffusionTensor3D< double >, 3 > TensorImageType; + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + typedef itk::Matrix< double, 3, 3 > MatrixType; + + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + + bool SetShCoefficients(vnl_vector< double > shCoefficients, double b0); + vnl_matrix SetFiberDirection(GradientType& fiberDirection); + void SetFaRange(double min, double max){ m_FaRange.first = min; m_FaRange.second = max; } + void SetAdcRange(double min, double max){ m_AdcRange.first = min; m_AdcRange.second = max; } + void SetMaxNumKernels(unsigned int max){ m_MaxNumKernels = max; } + unsigned int GetNumberOfKernels(); + std::pair< double, double > GetFaRange(){ return m_FaRange; } + std::pair< double, double > GetAdcRange(){ return m_AdcRange; } + unsigned int GetMaxNumKernels(){ return m_MaxNumKernels; } + void Clear(); + + std::vector< vnl_vector< double > > GetShCoefficients(){ return m_ShCoefficients; } + std::vector< double > GetB0Signals(){ return m_B0Signal; } + unsigned int GetShOrder(){ return m_ShOrder; } + int GetModelIndex(){ return m_ModelIndex; } + + double GetBaselineSignal(int index){ return m_B0Signal.at(index); } + vnl_vector< double > GetCoefficients(int listIndex){ return m_ShCoefficients.at(listIndex); } + + bool SampleKernels(Image::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage=nullptr, QballFilterType::CoefficientImageType::Pointer coeffImage=nullptr, ItkDoubleImageType::Pointer adcImage=nullptr); protected: - vnl_matrix Cart2Sph( GradientListType& gradients ); - void RandomModel(); - - std::vector< vnl_vector< double > > m_ShCoefficients; - std::vector< double > m_B0Signal; - std::vector< GradientType > m_PrototypeMaxDirection; - std::pair< double, double > m_AdcRange; - std::pair< double, double > m_FaRange; - unsigned int m_ShOrder; - int m_ModelIndex; - unsigned int m_MaxNumKernels; + vnl_matrix Cart2Sph( GradientListType& gradients ); + void RandomModel(); + + std::vector< vnl_vector< double > > m_ShCoefficients; + std::vector< double > m_B0Signal; + std::vector< GradientType > m_PrototypeMaxDirection; + std::pair< double, double > m_AdcRange; + std::pair< double, double > m_FaRange; + unsigned int m_ShOrder; + int m_ModelIndex; + unsigned int m_MaxNumKernels; }; } #include "mitkRawShModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h index 8e2ab01f8d..23e68f666e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h @@ -1,69 +1,70 @@ /*=================================================================== 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 _MITK_StickModel_H #define _MITK_StickModel_H #include namespace mitk { /** * \brief Generates the diffusion signal using an idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class StickModel : public DiffusionSignalModel< ScalarType > { public: - StickModel(); - template< class OtherType >StickModel(StickModel* model) - { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); - - this->m_BValue = model->GetBvalue(); - this->m_Diffusivity = model->GetDiffusivity(); - } - ~StickModel(); - - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); - - void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant - double GetDiffusivity() { return m_Diffusivity; } + StickModel(); + template< class OtherType >StickModel(StickModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); + + this->m_BValue = model->GetBvalue(); + this->m_Diffusivity = model->GetDiffusivity(); + } + ~StickModel(); + + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + + void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant + double GetDiffusivity() { return m_Diffusivity; } protected: - double m_Diffusivity; ///< Scalar diffusion constant + double m_Diffusivity; ///< Scalar diffusion constant }; } #include "mitkStickModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h index 3d8deb420d..f9522ba450 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h @@ -1,82 +1,83 @@ /*=================================================================== 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 _MITK_TensorModel_H #define _MITK_TensorModel_H #include #include namespace mitk { /** * \brief Generates diffusion measurement employing a second rank tensor model: e^(-bg^TDg) * */ template< class ScalarType = double > class TensorModel : public DiffusionSignalModel< ScalarType > { public: - TensorModel(); - template< class OtherType >TensorModel(TensorModel* model) - { - this->m_CompartmentId = model->m_CompartmentId; - this->m_T2 = model->GetT2(); - this->m_GradientList = model->GetGradientList(); - this->m_VolumeFractionImage = model->GetVolumeFractionImage(); - this->m_RandGen = model->GetRandomGenerator(); - - this->m_BValue = model->GetBvalue(); - this->m_KernelDirection = model->GetKernelDirection(); - this->m_KernelTensorMatrix = model->GetKernelTensorMatrix(); - } - ~TensorModel(); - - typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; - typedef itk::DiffusionTensor3D< ScalarType > ItkTensorType; - typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; - typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; - - /** Actual signal generation **/ - PixelType SimulateMeasurement(GradientType& fiberDirection); - ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); - - void SetDiffusivity1(double d1){ m_KernelTensorMatrix[0][0] = d1; } - void SetDiffusivity2(double d2){ m_KernelTensorMatrix[1][1] = d2; } - void SetDiffusivity3(double d3){ m_KernelTensorMatrix[2][2] = d3; } - double GetDiffusivity1() { return m_KernelTensorMatrix[0][0]; } - double GetDiffusivity2() { return m_KernelTensorMatrix[1][1]; } - double GetDiffusivity3() { return m_KernelTensorMatrix[2][2]; } - - GradientType GetKernelDirection(){ return m_KernelDirection; } - vnl_matrix_fixed GetKernelTensorMatrix(){ return m_KernelTensorMatrix; } + TensorModel(); + template< class OtherType >TensorModel(TensorModel* model) + { + this->m_CompartmentId = model->m_CompartmentId; + this->m_T1 = model->GetT1(); + this->m_T2 = model->GetT2(); + this->m_GradientList = model->GetGradientList(); + this->m_VolumeFractionImage = model->GetVolumeFractionImage(); + this->m_RandGen = model->GetRandomGenerator(); + + this->m_BValue = model->GetBvalue(); + this->m_KernelDirection = model->GetKernelDirection(); + this->m_KernelTensorMatrix = model->GetKernelTensorMatrix(); + } + ~TensorModel(); + + typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; + typedef itk::DiffusionTensor3D< ScalarType > ItkTensorType; + typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; + typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; + + /** Actual signal generation **/ + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); + + void SetDiffusivity1(double d1){ m_KernelTensorMatrix[0][0] = d1; } + void SetDiffusivity2(double d2){ m_KernelTensorMatrix[1][1] = d2; } + void SetDiffusivity3(double d3){ m_KernelTensorMatrix[2][2] = d3; } + double GetDiffusivity1() { return m_KernelTensorMatrix[0][0]; } + double GetDiffusivity2() { return m_KernelTensorMatrix[1][1]; } + double GetDiffusivity3() { return m_KernelTensorMatrix[2][2]; } + + GradientType GetKernelDirection(){ return m_KernelDirection; } + vnl_matrix_fixed GetKernelTensorMatrix(){ return m_KernelTensorMatrix; } protected: - /** Calculates tensor matrix from FA and ADC **/ - void UpdateKernelTensor(); - GradientType m_KernelDirection; ///< Direction of the kernel tensors principal eigenvector - vnl_matrix_fixed m_KernelTensorMatrix; ///< 3x3 matrix containing the kernel tensor values + /** Calculates tensor matrix from FA and ADC **/ + void UpdateKernelTensor(); + GradientType m_KernelDirection; ///< Direction of the kernel tensors principal eigenvector + vnl_matrix_fixed m_KernelTensorMatrix; ///< 3x3 matrix containing the kernel tensor values }; } #include "mitkTensorModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp index d0b152e50c..b75d4fcc51 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp @@ -1,333 +1,362 @@ /*=================================================================== 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. ===================================================================*/ #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include using namespace mitk; float CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { typedef itk::VectorImage< short, 3 > DwiImageType; try{ itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); unsigned int count = 0; float difference = 0; while(!it1.IsAtEnd()) { for (unsigned int i=0; iGetVectorLength(); ++i) { difference += abs(it1.Get()[i]-it2.Get()[i]); count++; } ++it1; ++it2; } return difference/count; } catch(...) { return -1; } return -1; } -FiberfoxParameters MakeProposal(FiberfoxParameters old_params, double temperature) +FiberfoxParameters MakeProposalRelaxation(FiberfoxParameters old_params, double temperature) { std::random_device r; std::default_random_engine randgen(r()); - std::uniform_int_distribution uint1(0, 2); + std::uniform_int_distribution uint1(0, 4); FiberfoxParameters new_params(old_params); int prop = uint1(randgen); -// prop = 2; switch(prop) { case 0: { std::normal_distribution normal_dist(0, new_params.m_SignalGen.m_SignalScale*0.1*temperature); float add = 0; while (add == 0) add = normal_dist(randgen); new_params.m_SignalGen.m_SignalScale += add; MITK_INFO << "Proposal Signal Scale: " << add << " (" << new_params.m_SignalGen.m_SignalScale << ")"; break; } case 1: { int model_index = rand()%new_params.m_NonFiberModelList.size(); double t2 = new_params.m_NonFiberModelList[model_index]->GetT2(); std::normal_distribution normal_dist(0, t2*0.1*temperature); double add = 0; while (add == 0) add = normal_dist(randgen); t2 += add; new_params.m_NonFiberModelList[model_index]->SetT2(t2); MITK_INFO << "Proposal T2 (Non-Fiber " << model_index << "): " << add << " (" << t2 << ")"; break; } case 2: { int model_index = rand()%new_params.m_FiberModelList.size(); double t2 = new_params.m_FiberModelList[model_index]->GetT2(); std::normal_distribution normal_dist(0, t2*0.1*temperature); double add = 0; while (add == 0) add = normal_dist(randgen); t2 += add; new_params.m_FiberModelList[model_index]->SetT2(t2); MITK_INFO << "Proposal T2 (Fiber " << model_index << "): " << add << " (" << t2 << ")"; break; } + case 3: + { + int model_index = rand()%new_params.m_NonFiberModelList.size(); + double t1 = new_params.m_NonFiberModelList[model_index]->GetT1(); + MITK_INFO << "T1: " << t1; + std::normal_distribution normal_dist(0, t1*0.1*temperature); + + double add = 0; + while (add == 0) + add = normal_dist(randgen); + t1 += add; + new_params.m_NonFiberModelList[model_index]->SetT1(t1); + MITK_INFO << "Proposal T1 (Non-Fiber " << model_index << "): " << add << " (" << t1 << ")"; + break; + } + case 4: + { + int model_index = rand()%new_params.m_FiberModelList.size(); + double t1 = new_params.m_FiberModelList[model_index]->GetT1(); + MITK_INFO << "T1: " << t1; + std::normal_distribution normal_dist(0, t1*0.1*temperature); + + double add = 0; + while (add == 0) + add = normal_dist(randgen); + t1 += add; + new_params.m_FiberModelList[model_index]->SetT1(t1); + MITK_INFO << "Proposal T1 (Fiber " << model_index << "): " << add << " (" << t1 << ")"; + break; + } } return new_params; } double UpdateDiffusivity(double d, double temperature) { std::random_device r; std::default_random_engine randgen(r()); std::normal_distribution normal_dist(0, d*0.1*temperature); double add = 0; while (add == 0) add = normal_dist(randgen); d += add; return d; } void ProposeDiffusivities(mitk::DiffusionSignalModel<>* signalModel, double temperature) { if (dynamic_cast*>(signalModel)) { mitk::StickModel<>* m = dynamic_cast*>(signalModel); m->SetDiffusivity(UpdateDiffusivity(m->GetDiffusivity(), temperature)); MITK_INFO << "d: " << m->GetDiffusivity(); } else if (dynamic_cast*>(signalModel)) { mitk::TensorModel<>* m = dynamic_cast*>(signalModel); m->SetDiffusivity1(UpdateDiffusivity(m->GetDiffusivity1(), temperature)); double new_d2 = UpdateDiffusivity(m->GetDiffusivity2(), temperature); m->SetDiffusivity2(new_d2); m->SetDiffusivity3(new_d2); MITK_INFO << "d1: " << m->GetDiffusivity1(); MITK_INFO << "d2: " << m->GetDiffusivity2(); MITK_INFO << "d3: " << m->GetDiffusivity3(); } else if (dynamic_cast*>(signalModel)) { mitk::BallModel<>* m = dynamic_cast*>(signalModel); m->SetDiffusivity(UpdateDiffusivity(m->GetDiffusivity(), temperature)); MITK_INFO << "d: " << m->GetDiffusivity(); } else if (dynamic_cast*>(signalModel)) { mitk::AstroStickModel<>* m = dynamic_cast*>(signalModel); m->SetDiffusivity(UpdateDiffusivity(m->GetDiffusivity(), temperature)); MITK_INFO << "d: " << m->GetDiffusivity(); } } FiberfoxParameters MakeProposalDiff(FiberfoxParameters old_params, double temperature) { std::random_device r; std::default_random_engine randgen(r()); std::uniform_int_distribution uint1(0, 1); FiberfoxParameters new_params(old_params); int prop = uint1(randgen); switch(prop) { case 0: { int model_index = rand()%new_params.m_NonFiberModelList.size(); ProposeDiffusivities( new_params.m_NonFiberModelList[model_index], temperature ); MITK_INFO << "Proposal D (Non-Fiber " << model_index << ")"; break; } case 1: { int model_index = rand()%new_params.m_FiberModelList.size(); ProposeDiffusivities( new_params.m_FiberModelList[model_index], temperature ); MITK_INFO << "Proposal D (Fiber " << model_index << ")"; break; } } return new_params; } /*! * \brief Command line interface to Fiberfox. * Simulate a diffusion-weighted image from a tractogram using the specified parameter file. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("FiberfoxOptimization"); parser.setCategory("Optimize Fiberfox Parameters"); parser.setContributor("MIC"); parser.setDescription("Command line interface to Fiberfox." " Simulate a diffusion-weighted image from a tractogram using the specified parameter file."); parser.setArgumentPrefix("--", "-"); parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameter file:", "fiberfox parameter file (.ffp)", us::Any(), false); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "Input tractogram or diffusion-weighted image.", us::Any(), false); parser.addArgument("template", "t", mitkCommandLineParser::String, "Template image:", "Use parameters of the template diffusion-weighted image.", us::Any(), false); parser.addArgument("iterations", "", mitkCommandLineParser::Int, "Iterations:", "Number of optimizations steps", 1000); parser.addArgument("start_temp", "", mitkCommandLineParser::Float, "Start temperature:", "", 1.0); parser.addArgument("end_temp", "", mitkCommandLineParser::Float, "End temperature:", "", 0.1); parser.addArgument("optimize_diff", "", mitkCommandLineParser::Bool, "Optimize diffusivities:", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string paramName = us::any_cast(parsedArgs["parameters"]); std::string input = us::any_cast(parsedArgs["input"]); int iterations=1000; if (parsedArgs.count("iterations")) iterations = us::any_cast(parsedArgs["iterations"]); float start_temp=1.0; if (parsedArgs.count("start_temp")) start_temp = us::any_cast(parsedArgs["start_temp"]); float end_temp=0.1; if (parsedArgs.count("end_temp")) end_temp = us::any_cast(parsedArgs["end_temp"]); bool optimize_diff=false; if (parsedArgs.count("optimize_diff")) optimize_diff = true; FiberfoxParameters parameters; parameters.LoadParameters(paramName); MITK_INFO << "Loading template image"; typedef itk::VectorImage< short, 3 > ItkDwiType; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {}); mitk::Image::Pointer dwi = dynamic_cast(mitk::IOUtil::Load(us::any_cast(parsedArgs["template"]), &functor)[0].GetPointer()); ItkDwiType::Pointer reference = mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi); parameters.m_SignalGen.m_ImageRegion = reference->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = reference->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = reference->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = reference->GetDirection(); parameters.SetBvalue(static_cast(dwi->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); parameters.SetGradienDirections(static_cast( dwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input, &functor)[0]; itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); tractsToDwiFilter->SetFiberBundle(dynamic_cast(inputData.GetPointer())); tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->Update(); ItkDwiType::Pointer sim = tractsToDwiFilter->GetOutput(); { mitk::Image::Pointer image = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.GetBvalue() ) ); mitk::DiffusionPropertyHelper propertyHelper( image ); propertyHelper.InitializeImage(); mitk::IOUtil::Save(image, "initial.dwi"); } MITK_INFO << "Iterations: " << iterations; MITK_INFO << "start_temp: " << start_temp; MITK_INFO << "end_temp: " << end_temp; double alpha = log(end_temp/start_temp); int accepted = 0; float last_diff = CompareDwi(sim, reference); for (int i=0; i<1000; ++i) { double temperature = start_temp * exp(alpha*(double)i/iterations); MITK_INFO << "Temperature: " << temperature << " (" << i+1 << "/" << iterations << ")"; FiberfoxParameters proposal; if (optimize_diff) proposal = MakeProposalDiff(parameters, temperature); else - proposal = MakeProposal(parameters, temperature); + proposal = MakeProposalRelaxation(parameters, temperature); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); tractsToDwiFilter->SetFiberBundle(dynamic_cast(inputData.GetPointer())); tractsToDwiFilter->SetParameters(proposal); tractsToDwiFilter->Update(); ItkDwiType::Pointer sim = tractsToDwiFilter->GetOutput(); std::cout.rdbuf (old); // <-- restore float diff = CompareDwi(sim, reference); if (last_diffGetOutput() ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.GetBvalue() ) ); mitk::DiffusionPropertyHelper propertyHelper( image ); propertyHelper.InitializeImage(); mitk::IOUtil::Save(image, "optimized.dwi"); - std::cout.rdbuf (old); // <-- restore proposal.SaveParameters("optimized.ffp"); + std::cout.rdbuf (old); // <-- restore accepted++; MITK_INFO << "Accepted proposal (acc. rate " << (float)accepted/(i+1) << ")"; parameters = FiberfoxParameters(proposal); last_diff = diff; } MITK_INFO << "\n\n\n"; } return EXIT_SUCCESS; }