ANIMA  4.0
animaODFResampleImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <animaODFFunctions.h>
5 
6 namespace anima
7 {
8 
9 template <typename TImageType, typename TInterpolatorPrecisionType>
10 void
13 {
14  if (!this->GetFiniteStrainReorientation())
15  itkExceptionMacro("ODF resampling only supports finite strain re-orientation");
16 
17  this->Superclass::BeforeThreadedGenerateData();
18 
19  unsigned int vectorSize = this->GetInput(0)->GetNumberOfComponentsPerPixel();
20  m_LOrder = (unsigned int)floor((-3 + sqrt(8*vectorSize + 1))/2);
21 
22  m_EulerAngles.resize(this->GetNumberOfWorkUnits());
23  m_ODFRotationMatrices.resize(this->GetNumberOfWorkUnits());
24 
25  for (unsigned int i = 0;i < this->GetNumberOfWorkUnits();++i)
26  m_EulerAngles[i].resize(3);
27 }
28 
29 template <typename TImageType, typename TInterpolatorPrecisionType>
30 void
32 ::ReorientInterpolatedModel(const InputPixelType &interpolatedModel, vnl_matrix <double> &modelOrientationMatrix,
33  InputPixelType &rotatedModel, itk::ThreadIdType threadId)
34 {
35  anima::GetEulerAnglesFromRotationMatrix(modelOrientationMatrix,m_EulerAngles[threadId]);
36  rotatedModel = interpolatedModel;
37 
38  for (unsigned int l = 2;l <= m_LOrder;l += 2)
39  {
40  anima::EstimateLocalODFRotationMatrix(m_ODFRotationMatrices[threadId],l,m_EulerAngles[threadId][0],
41  m_EulerAngles[threadId][1],m_EulerAngles[threadId][2]);
42 
43  unsigned int mBaseInd = (l*l + l + 2)/2 - l - 1;
44 
45  for (unsigned int m = 0;m <= 2*l;++m)
46  {
47  rotatedModel[mBaseInd + m] = 0;
48  for (unsigned int mp = 0;mp <= 2*l;++mp)
49  rotatedModel[mBaseInd + m] += m_ODFRotationMatrices[threadId](m,mp)*interpolatedModel[mBaseInd + mp];
50  }
51  }
52 }
53 
54 } // end namespace anima
virtual void BeforeThreadedGenerateData() ITK_OVERRIDE
Superclass::InputPixelType InputPixelType
virtual void ReorientInterpolatedModel(const InputPixelType &interpolatedModel, vnl_matrix< double > &modelOrientationMatrix, InputPixelType &rotatedModel, itk::ThreadIdType threadId) ITK_OVERRIDE
Needs to be implemented in sub-classes, does the actual re-orientation of the model.
void GetEulerAnglesFromRotationMatrix(vnl_matrix< double > &R, std::vector< double > &resVal)
void EstimateLocalODFRotationMatrix(vnl_matrix< double > &resVal, unsigned int l, double alpha, double beta, double gamma)