DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkSHCoefficientsRotation.h
Go to the documentation of this file.
1 
11 #ifndef __itkSHCoefficientsRotation_h
12 #define __itkSHCoefficientsRotation_h
13 
14 #include "itkObject.h"
15 #include "itkObjectFactory.h"
16 #include "utlNDArray.h"
17 #include "utlITKMacro.h"
18 #include "utl.h"
19 
20 namespace itk
21 {
22 
23 
34 template < class T = double >
35 class ITK_EXPORT SHCoefficientsRotation
36  : public Object
37 {
38 public:
41  typedef Object Superclass;
42  typedef SmartPointer<Self> Pointer;
43  typedef SmartPointer<const Self> ConstPointer;
44 
46  itkNewMacro(Self);
47 
49  itkTypeMacro( SHCoefficientsRotation, Object );
50 
51  typedef T PreciseType;
52 
55  typedef utl_shared_ptr<MatrixType> MatrixPointer;
56 
57  itkSetGetMacro(TessOrder, int);
58  itkSetGetMacro(MaxRank, int);
59 
60  void Initialize()
61  {
63  if (m_MaxRank/2>m_SHMatrixInverse->size())
64  {
66  MatrixPointer shMatrix = utl::ComputeSHMatrix(m_MaxRank, *m_Grad, CARTESIAN_TO_SPHERICAL);
67  m_SHMatrixInverse = utl_shared_ptr<std::vector<MatrixType> >(new std::vector<MatrixType>(m_MaxRank/2));
68  for ( int l = 2; l <= m_MaxRank; l += 2 )
69  {
70  MatrixType shMatrixL;
71  shMatrixL = shMatrix->GetNColumns(utl::RankToDimSH(l-2), 2*l+1);
72  (*m_SHMatrixInverse)[l/2-1] = utl::PInverseMatrix(shMatrixL, 1e-15);
73  }
74  }
75  }
76 
90  VectorType GetRotatedSHCoefficients(const VectorType& shInput, const MatrixType& rotationMatrix) const
91  {
92  utlSAException(m_SHMatrixInverse->size()==0).msg("need to use Initialize() first");
93  int rankReal = utl::DimToRankSH(shInput.Size());
94  utlSAException(rankReal>m_MaxRank)(rankReal)(m_MaxRank).msg("rank is larger than m_MaxRank.");
95  MatrixType gradt3Rotated = (*m_Grad) * rotationMatrix;
96  MatrixPointer shMatrixRotated = utl::ComputeSHMatrix(rankReal, gradt3Rotated,CARTESIAN_TO_SPHERICAL);
97  VectorType shRotated(shInput);
98  MatrixType tmp;
99 
100  for ( int l = 2; l <= rankReal; l += 2 )
101  {
102  int colstart = utl::RankToDimSH(l-2);
103  tmp = shMatrixRotated->GetNColumns(colstart, 2*l+1);
104  VectorType sfValueRotatedL = tmp* shInput.GetSubVector(utl::GetRange(colstart, colstart+2*l+1));
105  VectorType shRotatedL = (*m_SHMatrixInverse)[l/2-1]*sfValueRotatedL;
106  shRotated.SetSubVector(utl::GetRange(colstart, colstart+shRotatedL.Size()), shRotatedL);
107  }
108  return shRotated;
109  }
110 
111  static Pointer GetInstance()
112  {
113  if (!m_Instance)
114  m_Instance = Self::New();
115  return m_Instance;
116  }
117 
118 
119 protected:
120  SHCoefficientsRotation() : Superclass(),
121  m_Grad(new MatrixType()),
122  m_SHMatrixInverse(new std::vector<MatrixType>())
123  {
124  m_MaxRank = 10;
125  m_TessOrder = 3;
126  }
128 
129  void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE
130  {
131  Superclass::PrintSelf(os, indent);
132  PrintVar2(true, m_MaxRank, m_TessOrder, os <<indent);
133  utl::PrintUtlMatrix(*m_Grad, "m_Grad", " ", os << indent);
134  for ( int i = 0; i < m_SHMatrixInverse->size(); ++i )
135  {
136  PrintVar(true, os << indent, i);
137  utl::PrintUtlMatrix((*m_SHMatrixInverse)[i], "(*m_SHMatrixInverse)[i]", " ", os << indent);
138  }
139  }
140 
141  typename LightObject::Pointer InternalClone() const ITK_OVERRIDE
142  {
143  typename LightObject::Pointer loPtr = Superclass::InternalClone();
144  typename Self::Pointer rval = dynamic_cast<Self *>(loPtr.GetPointer());
145  if(rval.IsNull())
146  {
147  itkExceptionMacro(<< "downcast to type " << this->GetNameOfClass()<< " failed.");
148  }
149  rval->m_MaxRank = m_MaxRank;
150  rval->m_TessOrder = m_TessOrder;
151  rval->m_Grad = m_Grad;
152  rval->m_SHMatrixInverse = m_SHMatrixInverse;
153  return loPtr;
154  }
155 
156  static Pointer m_Instance;
157 
158  MatrixPointer m_Grad;
159  utl_shared_ptr<std::vector<MatrixType> > m_SHMatrixInverse;
160 
165 
166 private:
167  SHCoefficientsRotation(const Self&); //purposely not implemented
168  void operator=(const Self&); //purposely not implemented
169 
170 };
171 
172 template<>
174 
175 inline void
176 InitializeSHRotation(const int rank, const int tess)
177 {
178  typedef itk::SHCoefficientsRotation<double> SHRotationFilterType;
179  SHRotationFilterType::Pointer shRotate = SHRotationFilterType::GetInstance();
180  shRotate->SetTessOrder(tess);
181  shRotate->SetMaxRank(rank);
182  shRotate->Initialize();
183 }
184 
185 }
186 
187 
188 #endif
NDArray is a N-Dimensional array class (row-major, c version)
Definition: utlFunctors.h:131
static MatrixPointer GetGrad(const int tessorder, const int NoSymmetricDuple=DIRECTION_NODUPLICATE, const int mode=CARTESIAN_TO_SPHERICAL, const int flipx=DIRECTION_NOFLIP, const int flipy=DIRECTION_NOFLIP, const int flipz=DIRECTION_NOFLIP, const bool need_normalize=true)
helper functions specifically used in dmritool
void PrintUtlMatrix(const NDArray< T, 2 > &mat, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout)
STL namespace.
void InitializeSHRotation(const int rank, const int tess)
int DimToRankSH(const int dimm)
Definition: utlDMRI.h:182
utl::NDArray< double, 1 > VectorType
utl::NDArray< double, 2 > MatrixType
#define ITK_OVERRIDE
Definition: utlITKMacro.h:46
std::shared_ptr< NDArray< T, 2 > > ComputeSHMatrix(const unsigned int rank, const NDArray< T, 2 > &grad, const int mode)
Definition: utl.h:171
void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
SmartPointer< const Self > ConstPointer
LightObject::Pointer InternalClone() const ITK_OVERRIDE
#define PrintVar(cond, os,...)
Definition: utlCoreMacro.h:242
NDArray< T, 2 > PInverseMatrix(const NDArray< T, 2 > &mat, const double eps=1e-10)
#define utlSAException(expr)
Definition: utlCoreMacro.h:543
VectorType GetRotatedSHCoefficients(const VectorType &shInput, const MatrixType &rotationMatrix) const
GetRotatedSHCoefficient rotate SH coefficient vector by a given rotation matrix.
rotate SH coefficient vector by a given rotation matrix.
int RankToDimSH(const int shRank)
Definition: utlDMRI.h:193
#define itkSetGetMacro(name, type)
Definition: utlITKMacro.h:134
static void Initialize(const int tessorder)
utl_shared_ptr< MatrixType > MatrixPointer
utl_shared_ptr< std::vector< MatrixType > > m_SHMatrixInverse
SizeType Size() const
Definition: utlNDArray.h:321
std::vector< int > GetRange(const int start, const int end, const int space=1)
Definition: utlCore.h:1642
#define PrintVar2(cond, var1, var2, os)
Definition: utlCoreMacro.h:454