ANIMA  4.0
animaMCMEstimatorImageFilter.h
Go to the documentation of this file.
1 #pragma once
2 #include <cmath>
3 #include <random>
4 
6 #include <animaMCMImage.h>
7 #include <itkImage.h>
8 #include <itkSingleValuedCostFunction.h>
9 
11 #include <itkCostFunction.h>
12 #include <itkNonLinearOptimizer.h>
13 
15 #include <animaMCMConstants.h>
16 
17 namespace anima
18 {
19 
20 template <class InputPixelType, class OutputPixelType>
22  public anima::MaskedImageToImageFilter< itk::Image<InputPixelType,3>, anima::MCMImage<OutputPixelType,3> >
23 {
24 public:
27  typedef itk::Image<InputPixelType,3> TInputImage;
29  typedef itk::Image<InputPixelType,4> Image4DType;
31  typedef itk::Image<OutputPixelType,3> OutputScalarImageType;
32  typedef typename OutputScalarImageType::Pointer OutputScalarImagePointer;
34  typedef itk::SmartPointer<Self> Pointer;
35  typedef itk::SmartPointer<const Self> ConstPointer;
36 
37  typedef itk::Image<unsigned char,3> MoseImageType;
38  typedef typename MoseImageType::Pointer MoseImagePointer;
39 
40  typedef itk::VectorImage<OutputPixelType,3> VectorImageType;
41  typedef typename VectorImageType::Pointer VectorImagePointer;
42 
49 
50  typedef itk::NonLinearOptimizer OptimizerType;
51  typedef OptimizerType::Pointer OptimizerPointer;
52  typedef OptimizerType::ParametersType ParametersType;
53  typedef itk::CostFunction CostFunctionBaseType;
54  typedef CostFunctionBaseType::Pointer CostFunctionBasePointer;
55 
58  {
59  Gaussian = 0,
61  };
62 
64  {
65  Marginal = 0,
68  };
69 
71  itkNewMacro(Self)
72 
73 
75 
76 
77  typedef TInputImage InputImageType;
78  typedef TOutputImage OutputImageType;
79  typedef typename InputImageType::Pointer InputImagePointer;
80  typedef typename OutputImageType::PixelType VariableLengthVectorType;
81  typedef typename OutputImageType::Pointer OutputImagePointer;
82 
84  typedef typename Superclass::MaskImageType MaskImageType;
85  typedef typename Superclass::InputImageRegionType InputImageRegionType;
86  typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
87 
88  typedef vnl_vector_fixed <double,3> GradientType;
89 
90  // Acquisition-related parameters
91  void SetGradientStrengths(std::vector <double> &mb) {m_GradientStrengths = mb;}
92  itkSetMacro(SmallDelta, double)
93  itkSetMacro(BigDelta, double)
94  void AddGradientDirection(unsigned int i, GradientType &grad);
95  itkSetMacro(B0Threshold, double)
96 
97  // Optimizer-related parameters
98  void SetOptimizer(std::string &opt) {m_Optimizer = opt;}
99  itkSetMacro(AbsoluteCostChange, double)
100 
101  itkSetMacro(MLEstimationStrategy, MaximumLikelihoodEstimationMode)
102  itkGetMacro(MLEstimationStrategy, MaximumLikelihoodEstimationMode)
103 
104  // Model-related parameters
105  itkSetMacro(ModelWithFreeWaterComponent, bool)
106  itkSetMacro(ModelWithStationaryWaterComponent, bool)
107  itkSetMacro(ModelWithRestrictedWaterComponent, bool)
108  itkSetMacro(ModelWithStaniszComponent, bool)
109 
110  itkSetMacro(NoiseType, SignalNoiseType)
111  itkGetMacro(NoiseType, SignalNoiseType)
112  itkSetMacro(CompartmentType, CompartmentType)
113 
114  itkSetMacro(NumberOfCoils, unsigned int)
115  itkGetMacro(NumberOfCoils, unsigned int)
116  itkSetMacro(NumberOfCompartments, unsigned int)
117  itkSetMacro(FindOptimalNumberOfCompartments, bool)
118 
119  itkSetMacro(UseConstrainedDiffusivity, bool)
120  itkSetMacro(UseConstrainedFreeWaterDiffusivity, bool)
121  itkSetMacro(UseConstrainedIRWDiffusivity, bool)
122  itkSetMacro(UseConstrainedStaniszDiffusivity, bool)
123  itkSetMacro(UseConstrainedStaniszRadius, bool)
124 
125  itkSetMacro(UseConstrainedOrientationConcentration, bool)
126  itkSetMacro(UseConstrainedExtraAxonalFraction, bool)
127  itkSetMacro(UseCommonConcentrations, bool)
128  itkSetMacro(UseCommonExtraAxonalFractions, bool)
129 
130  itkSetMacro(UseCommonDiffusivities, bool)
131 
132  std::string GetOptimizer() {return m_Optimizer;}
133 
134  std::vector <double> & GetGradientStrengths() {return m_GradientStrengths;}
135  itkGetMacro(SmallDelta, double)
136  itkGetMacro(BigDelta, double)
137  std::vector< GradientType > &GetGradientDirections() {return m_GradientDirections;}
138 
139  MCMCreatorType *GetMCMCreator(unsigned int i) {return m_MCMCreators[i];}
141 
142  // Output options
143  OutputScalarImageType *GetAICcVolume () {return m_AICcVolume;}
144  OutputScalarImageType *GetB0Volume () {return m_B0Volume;}
145  OutputScalarImageType *GetSigmaSquareVolume () {return m_SigmaSquareVolume;}
146 
147  void SetMoseVolume (MoseImageType *img) {m_MoseVolume = img; m_ExternalMoseVolume = true;}
148  MoseImageType *GetMoseVolume () {return m_MoseVolume;}
149 
150  void WriteMCMOutput(std::string fileName);
151 
152  itkSetMacro(AxialDiffusivityValue, double)
153  itkSetMacro(StaniszDiffusivityValue, double)
154  itkSetMacro(IRWDiffusivityValue, double)
155  itkSetMacro(RadialDiffusivity1Value, double)
156  itkSetMacro(RadialDiffusivity2Value, double)
157 
158  itkSetMacro(XTolerance, double)
159  itkSetMacro(FTolerance, double)
160  itkSetMacro(MaxEval, unsigned int)
161 
162 protected:
164  {
165  m_AICcVolume = 0;
166  m_B0Volume = 0;
167  m_SigmaSquareVolume = 0;
168  m_MoseVolume = 0;
169 
170  m_GradientStrengths.clear();
171  m_GradientDirections.clear();
172 
173  m_NumberOfDictionaryEntries = 500;
174  m_Optimizer = "bobyqa";
175  m_AbsoluteCostChange = 0.01;
176  m_B0Threshold = 0;
177  m_MLEstimationStrategy = Marginal;
178 
179  m_ModelWithFreeWaterComponent = true;
180  m_ModelWithStationaryWaterComponent = true;
181  m_ModelWithRestrictedWaterComponent = true;
182  m_ModelWithStaniszComponent = true;
183 
184  m_NoiseType = Gaussian;
185  m_CompartmentType = anima::Tensor;
186 
187  m_NumberOfCoils = 1;
188  m_NumberOfCompartments = 2;
189  m_FindOptimalNumberOfCompartments = true;
190 
191  m_UseConstrainedDiffusivity = false;
192  m_UseConstrainedFreeWaterDiffusivity = true;
193  m_UseConstrainedIRWDiffusivity = true;
194  m_UseConstrainedStaniszDiffusivity = true;
195  m_UseConstrainedStaniszRadius = true;
196  m_UseCommonDiffusivities = false;
197 
198  m_UseConstrainedOrientationConcentration = false;
199  m_UseConstrainedExtraAxonalFraction = false;
200  m_UseCommonConcentrations = false;
201  m_UseCommonExtraAxonalFractions = false;
202 
203  m_AxialDiffusivityValue = 1.71e-3;
204  m_StaniszDiffusivityValue = 1.71e-3;
205  m_IRWDiffusivityValue = 7.5e-4;
206  m_RadialDiffusivity1Value = 1.9e-4;
207  m_RadialDiffusivity2Value = 1.5e-4;
208 
209  m_NumberOfImages = 0;
210  m_ExternalMoseVolume = false;
211 
212  m_MaxEval = 0;
213  m_XTolerance = 0;
214  m_FTolerance = 0;
215 
216  m_SmallDelta = anima::DiffusionSmallDelta;
217  m_BigDelta = anima::DiffusionBigDelta;
218  }
219 
221  {
222  for (unsigned int i = 0;i < m_MCMCreators.size();++i)
223  delete m_MCMCreators[i];
224  m_MCMCreators.clear();
225  }
226 
227  void CheckComputationMask() ITK_OVERRIDE;
228 
229  void GenerateOutputInformation() ITK_OVERRIDE;
230  virtual void BeforeThreadedGenerateData() ITK_OVERRIDE;
231  void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE;
232 
234  virtual CostFunctionBasePointer CreateCostFunction(std::vector<double> &observedSignals, MCMPointer &mcmModel);
235 
237  OptimizerPointer CreateOptimizer(CostFunctionBasePointer &cost, itk::Array<double> &lowerBounds, itk::Array<double> &upperBounds);
238 
240  void EstimateFreeWaterModel(MCMPointer &mcmValue, std::vector <double> &observedSignals, itk::ThreadIdType threadId,
241  double &aiccValue, double &b0Value, double &sigmaSqValue);
242 
244  void OptimizeNonIsotropicCompartments(MCMPointer &mcmValue, unsigned int currentNumberOfCompartments,
245  std::vector <double> &observedSignals, itk::ThreadIdType threadId,
246  double &aiccValue, double &b0Value, double &sigmaSqValue);
247 
249  void InitialOrientationsEstimation(MCMPointer &mcmValue, bool authorizedNegativeB0Value, unsigned int currentNumberOfCompartments,
250  std::vector <double> &observedSignals, itk::ThreadIdType threadId,
251  double &aiccValue, double &b0Value, double &sigmaSqValue);
252 
254  void ModelEstimation(MCMPointer &mcmValue, bool authorizedNegativeB0Value, std::vector <double> &observedSignals,
255  itk::ThreadIdType threadId, double &aiccValue, double &b0Value, double &sigmaSqValue);
256 
258  double PerformSingleOptimization(ParametersType &p, CostFunctionBasePointer &cost, itk::Array<double> &lowerBounds,
259  itk::Array<double> &upperBounds);
260 
262  virtual void SparseInitializeSticks(MCMPointer &complexModel, bool authorizeNegativeB0Value,
263  std::vector<double> &observedSignals, itk::ThreadIdType threadId);
264 
266  virtual void InitializeModelFromSimplifiedOne(MCMPointer &simplifiedModel, MCMPointer &complexModel);
267 
270 
272  void GetProfiledInformation(CostFunctionBasePointer &cost, MCMPointer &mcm, double &b0Value, double &sigmaSqValue);
273 
275  double ComputeAICcValue(MCMPointer &mcmValue, double costValue);
276 
279 
282 
285  MCMType::ListType &workVec,ParametersType &p);
286 
289  MCMType::ListType &workVec,ParametersType &p);
290 
291 private:
292  ITK_DISALLOW_COPY_AND_ASSIGN(MCMEstimatorImageFilter);
293 
295  void InitializeDictionary();
296 
297  double m_SmallDelta, m_BigDelta;
298  std::vector <double> m_GradientStrengths;
299  std::vector< GradientType > m_GradientDirections;
300 
301  OutputScalarImagePointer m_B0Volume;
302  OutputScalarImagePointer m_SigmaSquareVolume;
303  OutputScalarImagePointer m_AICcVolume;
304  MoseImagePointer m_MoseVolume;
305 
306  std::vector <MCMCreatorType *> m_MCMCreators;
307 
308  std::string m_Optimizer;
309 
311  vnl_matrix <double> m_SparseSticksDictionary;
312  unsigned int m_NumberOfDictionaryEntries;
313  std::vector < std::vector <double> > m_DictionaryDirections;
314 
315  double m_B0Threshold;
316  unsigned int m_NumberOfImages;
317 
318  double m_AbsoluteCostChange;
319  MaximumLikelihoodEstimationMode m_MLEstimationStrategy;
320 
321  bool m_ModelWithFreeWaterComponent, m_ModelWithStationaryWaterComponent, m_ModelWithRestrictedWaterComponent, m_ModelWithStaniszComponent;
322 
323  SignalNoiseType m_NoiseType;
324  CompartmentType m_CompartmentType;
325 
326  unsigned int m_NumberOfCoils;
327  unsigned int m_NumberOfCompartments;
328  bool m_FindOptimalNumberOfCompartments;
329 
330  bool m_UseConstrainedDiffusivity;
331  bool m_UseConstrainedFreeWaterDiffusivity;
332  bool m_UseConstrainedIRWDiffusivity;
333  bool m_UseConstrainedStaniszRadius;
334  bool m_UseConstrainedStaniszDiffusivity;
335  bool m_UseCommonDiffusivities;
336 
337  bool m_UseConstrainedOrientationConcentration;
338  bool m_UseConstrainedExtraAxonalFraction;
339  bool m_UseCommonConcentrations;
340  bool m_UseCommonExtraAxonalFractions;
341 
342  double m_AxialDiffusivityValue;
343  double m_IRWDiffusivityValue;
344  double m_StaniszDiffusivityValue;
345  double m_RadialDiffusivity1Value;
346  double m_RadialDiffusivity2Value;
347 
348  bool m_ExternalMoseVolume;
349 
350  unsigned int m_MaxEval;
351  double m_XTolerance;
352  double m_FTolerance;
353 
355  std::vector < std::vector <double> > m_ValuesCoarseGrid;
356 };
357 
358 } // end namespace anima
359 
void TensorCoarseGridInitialization(MCMPointer &mcmUpdateValue, CostFunctionBasePointer &cost, MCMType::ListType &workVec, ParametersType &p)
Coarse grid initialization of tensor model.
OptimizerType::ParametersType ParametersType
anima::MultiCompartmentModelCreator::CompartmentType CompartmentType
std::vector< GradientType > & GetGradientDirections()
std::vector< double > & GetGradientStrengths()
virtual void InitializeModelFromSimplifiedOne(MCMPointer &simplifiedModel, MCMPointer &complexModel)
Performs initialization from simplified model with the same number of compartments.
anima::MultiCompartmentModelCreator MCMCreatorType
void EstimateFreeWaterModel(MCMPointer &mcmValue, std::vector< double > &observedSignals, itk::ThreadIdType threadId, double &aiccValue, double &b0Value, double &sigmaSqValue)
Specific method for N=0 compartments estimation (only free water)
double PerformSingleOptimization(ParametersType &p, CostFunctionBasePointer &cost, itk::Array< double > &lowerBounds, itk::Array< double > &upperBounds)
Performs an optimization of the supplied cost function and parameters using the specified optimizer(s...
itk::Image< unsigned char, 3 > MoseImageType
void SetGradientStrengths(std::vector< double > &mb)
STL namespace.
virtual void SparseInitializeSticks(MCMPointer &complexModel, bool authorizeNegativeB0Value, std::vector< double > &observedSignals, itk::ThreadIdType threadId)
Performs initialization from single DTI.
BaseCompartment::ListType ListType
SignalNoiseType
Denotes noise type on input signal. Based on this, a different cost function should be created...
MCMEstimatorImageFilter< InputPixelType, OutputPixelType > Self
void ComputeExtraAxonalAndKappaCoarseGrids()
Computes extra axonal and kappa coarse grids (used for NODDI and DDI initialization) ...
Superclass::InputImageRegionType InputImageRegionType
const double DiffusionBigDelta
Default big delta value (classical values)
void AddGradientDirection(unsigned int i, GradientType &grad)
const double DiffusionSmallDelta
Default small delta value (classical values)
CostFunctionBaseType::Pointer CostFunctionBasePointer
itk::Image< InputPixelType, 3 > TInputImage
double ComputeAICcValue(MCMPointer &mcmValue, double costValue)
Compute AICc value from a cost function value and model.
void GetProfiledInformation(CostFunctionBasePointer &cost, MCMPointer &mcm, double &b0Value, double &sigmaSqValue)
Utility function to get profiled data from the cost function into an MCM model.
itk::Image< OutputPixelType, 3 > OutputScalarImageType
DiffusionModelCompartmentType
void ModelEstimation(MCMPointer &mcmValue, bool authorizedNegativeB0Value, std::vector< double > &observedSignals, itk::ThreadIdType threadId, double &aiccValue, double &b0Value, double &sigmaSqValue)
Doing estimation, calling initialization procedure until ball and zeppelin, returns AICc value...
void OptimizeNonIsotropicCompartments(MCMPointer &mcmValue, unsigned int currentNumberOfCompartments, std::vector< double > &observedSignals, itk::ThreadIdType threadId, double &aiccValue, double &b0Value, double &sigmaSqValue)
Doing estimation of non isotropic compartments (for a given number of anisotropic compartments) ...
void ComputeTensorRadialDiffsAndAzimuthCoarseGrids()
Computes extra axonal and kappa coarse grids (used for tensor final initaialization) ...
anima::MCMImage< OutputPixelType, 3 > ModelImageType
itk::VectorImage< OutputPixelType, 3 > VectorImageType
itk::SmartPointer< const Self > ConstPointer
Superclass::OutputImageRegionType OutputImageRegionType
void ExtraAxonalAndKappaCoarseGridInitialization(MCMPointer &mcmUpdateValue, CostFunctionBasePointer &cost, MCMType::ListType &workVec, ParametersType &p)
Coarse grid initialization of NODDI and DDI models.
MultiCompartmentModel: holds several diffusion compartments, ordered by type It also handles weights ...
vnl_vector_fixed< double, 3 > GradientType
anima::MaskedImageToImageFilter< TInputImage, TOutputImage > Superclass
OptimizerPointer CreateOptimizer(CostFunctionBasePointer &cost, itk::Array< double > &lowerBounds, itk::Array< double > &upperBounds)
Create an optimizer following the optimizer type and estimation mode.
OutputImageType::PixelType VariableLengthVectorType
void InitialOrientationsEstimation(MCMPointer &mcmValue, bool authorizedNegativeB0Value, unsigned int currentNumberOfCompartments, std::vector< double > &observedSignals, itk::ThreadIdType threadId, double &aiccValue, double &b0Value, double &sigmaSqValue)
Doing estimation only of multiple orientations.
virtual MCMCreatorType * GetNewMCMCreatorInstance()
anima::MCMImage< OutputPixelType, 3 > TOutputImage
itk::Image< InputPixelType, 4 > Image4DType
double GetCostValue(CostFunctionBasePointer &cost, ParametersType &p)
Utility function to get a value from a cost function.
virtual void BeforeThreadedGenerateData() ITK_OVERRIDE
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
OutputScalarImageType::Pointer OutputScalarImagePointer
OutputScalarImageType * GetSigmaSquareVolume()
MCMCreatorType * GetMCMCreator(unsigned int i)
Really this class is some simplified factory that creates the MCM that it knows.
virtual CostFunctionBasePointer CreateCostFunction(std::vector< double > &observedSignals, MCMPointer &mcmModel)
Create a cost function following the noise type and estimation mode.