DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
utl4thOrderTensor.h
Go to the documentation of this file.
1 
10 #ifndef __utl4thOrderTensor_h
11 #define __utl4thOrderTensor_h
12 
13 #include "utlNDArray.h"
14 
15 namespace utl
16 {
17 
18 
19 template <class T>
20 inline void
22 
23 
32 template < class T >
33 class NDArray<T, 4> : public NDArrayBase<T,4>
34 {
35 public:
36  typedef NDArray Self;
38 
39  typedef typename Superclass::ValueType ValueType;
41 
42  typedef typename Superclass::SizeType SizeType;
43  typedef typename Superclass::ShapeType ShapeType;
44  typedef typename Superclass::Pointer Pointer;
46  typedef typename Superclass::Reference Reference;
48  typedef typename Superclass::Iterator Iterator;
50 
53 
54  using Superclass::SetData;
56  using Superclass::ReSize;
57  using Superclass::operator();
58  using Superclass::operator=;
59  using Superclass::operator+=;
60  using Superclass::operator-=;
61  using Superclass::operator%=;
62  using Superclass::operator/=;
63 
64  NDArray() : Superclass()
65  { }
66  NDArray(const NDArray<T,4>& mat) : Superclass(mat)
67  { }
68  template<typename EType>
69  NDArray(const Expr<EType, typename EType::ValueType>& expr) : Superclass(expr)
70  { }
71 
72  template< typename TMatrixValueType >
73  NDArray(const NDArray< TMatrixValueType,4> & r) : Superclass(r)
74  { }
75 
77  {
78  operator=(std::move(r));
79  }
80 
81  explicit NDArray(const ShapeType& shape) : Superclass(shape) { }
82 
83  NDArray(const T* vec, const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3) : Superclass()
84  {
85  SizeType shape[4];
86  shape[0]=s0, shape[1]=s1, shape[2]=s2, shape[3]=s3;
88  utl::cblas_copy<T>(this->Size(), vec, 1, this->m_Data, 1);
89  }
90 
91  NDArray(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3) : Superclass()
92  {
93  SizeType shape[4];
94  shape[0]=s0, shape[1]=s1, shape[2]=s2, shape[3]=s3;
96  }
97 
98  NDArray(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3, const T r) : Superclass()
99  {
100  SizeType shape[4];
101  shape[0]=s0, shape[1]=s1, shape[2]=s2, shape[3]=s3;
103  this->Fill(r);
104  }
105 
111  NDArray(const T* vec, const ShapeType& shape) : Superclass(vec,shape) { }
112 
116  NDArray(const ShapeType& shape, const T r) : Superclass(shape, r) { }
117 
119  {
121  return *this;
122  }
123 
125  {
126  if ( this != &r )
127  {
128  this->Clear();
129  this->Swap(r);
130  }
131  return *this;
132  }
133 
134  // 4th order tensor related
135 
136  UTL_ALWAYS_INLINE Reference operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
137  { return this->m_Data[i*this->m_OffSetTable[1] + j*this->m_OffSetTable[2] + k*this->m_OffSetTable[3]+l];}
138  UTL_ALWAYS_INLINE ConstReference operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
139  { return this->m_Data[i*this->m_OffSetTable[1] + j*this->m_OffSetTable[2] + k*this->m_OffSetTable[3]+l];}
140 
141  inline bool ReSize(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3)
142  {
143  SizeType shape[4];
144  shape[0]=s0, shape[1]=s1, shape[2]=s2, shape[3]=s3;
145  return Superclass::ReSize(shape);
146  }
147 
148  NDArray<T,2> GetRefSubMatrix(const int i, const int j) const
149  {
150  NDArray<T,2> mat;
151  unsigned index[4];
152  index[0]=i, index[1]=j, index[2]=0, index[3]=0;
153  unsigned shape[2];
154  shape[0]=this->m_Shape[2], shape[1]=this->m_Shape[3];
155  mat.SetData(this->m_Data+this->GetOffset(index), shape);
156  return mat;
157  }
158 
161  {
162  int size = this->m_Shape[0];
163  utlException(this->m_Shape[1]!=size || this->m_Shape[2]!=size || this->m_Shape[3]!=size, "should have the same size in every dimension");
164  SizeType index1[Dimension], index2[Dimension];
165 
166  NDArray<int, Dimension> isFilled(this->m_Shape, 0);
167  for ( int i = 0; i < size; ++i )
168  for ( int j = 0; j < size; ++j )
169  for ( int k = 0; k < size; ++k )
170  for ( int l = 0; l < size; ++l )
171  {
172  // utlPrintVar4(true, i,j,k,l);
173  index1[0]=i, index1[1]=j, index1[2]=k, index1[3]=l;
174  index2[0]=k, index2[1]=l, index2[2]=i, index2[3]=j;
175  if (isFilled(index1)==0 || isFilled(index2)==0)
176  {
177  T val = ((*this)(index1) + (*this)(index2))/2.0;
178  (*this)(index1)=val;
179  (*this)(index2)=val;
180 
181  isFilled(index1)=1;
182  isFilled(index2)=1;
183  }
184  }
185  }
186 
188  bool IsMajorSymmetric() const
189  {
190  int size = this->m_Shape[0];
191  utlException(this->m_Shape[1]!=size || this->m_Shape[2]!=size || this->m_Shape[3]!=size, "should have the same size in every dimension");
192  SizeType index1[Dimension], index2[Dimension];
193  NDArray<int, Dimension> isFilled(this->m_Shape, 0);
194  for ( int i = 0; i < size; ++i )
195  for ( int j = 0; j < size; ++j )
196  for ( int k = 0; k < size; ++k )
197  for ( int l = 0; l < size; ++l )
198  {
199  index1[0]=i, index1[1]=j, index1[2]=k, index1[3]=l;
200  index2[0]=k, index2[1]=l, index2[2]=i, index2[3]=j;
201  if (isFilled(index1)==0 || isFilled(index2)==0)
202  {
203  T val1 = (*this)(index1);
204  T val2 = (*this)(index2);
205  if (!utl::IsSame(val1, val2, 1e-10) )
206  {
207  // utl::PrintVector(index1, 4, "index1");
208  // utl::PrintVector(index2, 4, "index2");
209  // utlPrintVar2(true, (*this)(index1), (*this)(index2));
210  return false;
211  }
212  isFilled(index1)=1;
213  isFilled(index2)=1;
214  }
215  }
216  return true;
217  }
218 
220  {
221  int size = this->m_Shape[0];
222  utlException(this->m_Shape[1]!=size || this->m_Shape[2]!=size || this->m_Shape[3]!=size, "should have the same size in every dimension");
223  SizeType index1[Dimension], index2[Dimension], index3[Dimension], index4[Dimension];
224 
225  NDArray<int, Dimension> isFilled(this->m_Shape, 0);
226  for ( int i = 0; i < size; ++i )
227  for ( int j = 0; j < size; ++j )
228  for ( int k = 0; k < size; ++k )
229  for ( int l = 0; l < size; ++l )
230  {
231  // utlPrintVar4(true, i,j,k,l);
232  index1[0]=i, index1[1]=j, index1[2]=k, index1[3]=l;
233  index2[0]=j, index2[1]=i, index2[2]=k, index2[3]=l;
234  index3[0]=i, index3[1]=j, index3[2]=l, index3[3]=k;
235  index4[0]=j, index4[1]=i, index4[2]=l, index4[3]=k;
236 
237  if (isFilled(index1)==0 || isFilled(index2)==0 || isFilled(index3)==0 || isFilled(index4)==0)
238  {
239  T val1 = (*this)(index1),
240  val2 = (*this)(index2),
241  val3 = (*this)(index3),
242  val4 = (*this)(index4);
243  (*this)(index1) = (val1+val2+val3+val4)/4.0;
244  (*this)(index2) = (val1+val2+val3+val4)/4.0;
245  (*this)(index3) = (val1+val2+val3+val4)/4.0;
246  (*this)(index4) = (val1+val2+val3+val4)/4.0;
247 
248  isFilled(index1)=1;
249  isFilled(index2)=1;
250  isFilled(index3)=1;
251  isFilled(index4)=1;
252  }
253  }
254  }
255 
256  bool IsMinorSymmetric() const
257  {
258  int size = this->m_Shape[0];
259  utlException(this->m_Shape[1]!=size || this->m_Shape[2]!=size || this->m_Shape[3]!=size, "should have the same size in every dimension");
260 
261  NDArray<int, Dimension> isFilled(this->m_Shape, 0);
262  SizeType index1[Dimension], index2[Dimension], index3[Dimension], index4[Dimension];
263  for ( int i = 0; i < size; ++i )
264  for ( int j = 0; j < size; ++j )
265  for ( int k = 0; k < size; ++k )
266  for ( int l = 0; l < size; ++l )
267  {
268  index1[0]=i, index1[1]=j, index1[2]=k, index1[3]=l;
269  index2[0]=j, index2[1]=i, index2[2]=k, index2[3]=l;
270  index3[0]=i, index3[1]=j, index3[2]=l, index3[3]=k;
271  index4[0]=j, index4[1]=i, index4[2]=l, index4[3]=k;
272 
273  if (isFilled(index1)==0 || isFilled(index2)==0 || isFilled(index3)==0 || isFilled(index4)==0)
274  {
275  T val1 = (*this)(index1),
276  val2 = (*this)(index2),
277  val3 = (*this)(index3),
278  val4 = (*this)(index4);
279  if (!utl::IsSame(val1, val2, 1e-10) )
280  return false;
281  if (!utl::IsSame(val1, val3, 1e-10) )
282  return false;
283  if (!utl::IsSame(val1, val4, 1e-10) )
284  return false;
285  isFilled(index1)=1;
286  isFilled(index2)=1;
287  isFilled(index3)=1;
288  isFilled(index4)=1;
289  }
290  }
291  return true;
292  }
293 
295  {
296  utlException(this->m_Shape[0]!=3 || this->m_Shape[1]!=3 || this->m_Shape[2]!=3 || this->m_Shape[3]!=3, "should have size (3,3,3,3)");
297  NDArray<T, 2> mat;
298  utl::Convert4To2Tensor(*this, mat);
299  mat.EigenDecompositionNonSymmetricMatrix(valReal, valImg);
300  }
301 
302  void EigenDecompositionWithMinorSymmetry (NDArray<T,1>& valReal, NDArray<T,1>& valImg, std::vector<NDArray<T,2> >& matRealR, std::vector<NDArray<T,2> >& matImgR) const
303  {
304  utlException(this->m_Shape[0]!=3 || this->m_Shape[1]!=3 || this->m_Shape[2]!=3 || this->m_Shape[3]!=3, "should have size (3,3,3,3)");
305  matRealR.resize(6); matImgR.resize(6);
306  NDArray<T, 2> mat, vecRealR, vecImgR;
307  utl::Convert4To2Tensor(*this, mat);
308  mat.EigenDecompositionNonSymmetricMatrix(valReal, valImg, vecRealR, vecImgR);
309  NDArray<T,1> vecTmp(6);
310  for ( int i = 0; i < valReal.Size(); ++i )
311  {
312  for ( int j = 0; j < 6; ++j )
313  vecTmp[j] = vecRealR(j,i);
314  utl::Convert1To2Tensor(vecTmp, matRealR[i]);
315  for ( int j = 0; j < 6; ++j )
316  vecTmp[j] = vecImgR(j,i);
317  utl::Convert1To2Tensor(vecTmp, matImgR[i]);
318  }
319  }
320 
321  void EigenDecompositionWithMinorSymmetry (NDArray<T,1>& valReal, NDArray<T,1>& valImg, std::vector<NDArray<T,2> >& matRealR, std::vector<NDArray<T,2> >& matImgR, std::vector<NDArray<T,2> >& matRealL, std::vector<NDArray<T,2> >& matImgL) const
322  {
323  utlException(this->m_Shape[0]!=3 || this->m_Shape[1]!=3 || this->m_Shape[2]!=3 || this->m_Shape[3]!=3, "should have size (3,3,3,3)");
324  matRealR.resize(6); matImgR.resize(6);
325  matRealL.resize(6); matImgL.resize(6);
326  NDArray<T, 2> mat, vecRealR, vecImgR, vecRealL, vecImgL;
327  utl::Convert4To2Tensor(*this, mat);
328  mat.EigenDecompositionNonSymmetricMatrix(valReal, valImg, vecRealR, vecImgR, vecRealL, vecImgL);
329  NDArray<T,1> vecTmp(6);
330  for ( int i = 0; i < valReal.Size(); ++i )
331  {
332  for ( int j = 0; j < 6; ++j )
333  vecTmp[j] = vecRealR(j,i);
334  utl::Convert1To2Tensor(vecTmp, matRealR[i]);
335  for ( int j = 0; j < 6; ++j )
336  vecTmp[j] = vecImgR(j,i);
337  utl::Convert1To2Tensor(vecTmp, matImgR[i]);
338 
339  for ( int j = 0; j < 6; ++j )
340  vecTmp[j] = vecRealL(j,i);
341  utl::Convert1To2Tensor(vecTmp, matRealL[i]);
342  for ( int j = 0; j < 6; ++j )
343  vecTmp[j] = vecImgL(j,i);
344  utl::Convert1To2Tensor(vecTmp, matImgL[i]);
345  }
346  }
347 
348 };
349 
350 
351 }
352 
353 
354 #endif
const ValueType * ConstIterator
Definition: utlNDArray.h:171
Reference operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Superclass::Reference Reference
NDArray is a N-Dimensional array class (row-major, c version)
Definition: utlFunctors.h:131
NDArray(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3, const T r)
NDArray<T,1> is a vector class which uses blas mkl.
Definition: utlVector.h:36
void EigenDecompositionWithMinorSymmetry(NDArray< T, 1 > &valReal, NDArray< T, 1 > &valImg, std::vector< NDArray< T, 2 > > &matRealR, std::vector< NDArray< T, 2 > > &matImgR) const
base class for expression
Definition: utlExpression.h:36
SizeType m_Shape[Dimension]
Definition: utlNDArray.h:1168
bool ReSize(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3)
void EigenDecompositionWithMinorSymmetry(NDArray< T, 1 > &valReal, NDArray< T, 1 > &valImg) const
void Convert4To2Tensor(const utl::NDArray< T, 4 > &tensor, utl::NDArray< T, 2 > &mat)
NDArray(const T *vec, const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3)
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
static constexpr SizeType SubDimension
Definition: utlNDArray.h:162
Superclass::ConstPointer ConstPointer
NDArray<T,2> is a row-major matrix.
Definition: utlMatrix.h:37
bool ReSize(const ShapeType &shape)
Definition: utlNDArray.h:697
const ValueType & ConstReference
Definition: utlNDArray.h:179
Superclass::ConstReference ConstReference
ConstReference operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
NDArray(const Expr< EType, typename EType::ValueType > &expr)
void SetData(T *const data, const ShapeType &shape)
Definition: utlNDArray.h:677
NDArray< T, Dim > & operator=(const NDArray< T, Dim > &r)
Definition: utlNDArray.h:116
NDArray(const SizeType s0, const SizeType s1, const SizeType s2, const SizeType s3)
#define UTL_ALWAYS_INLINE
Definition: utlCoreMacro.h:76
NDArray(const ShapeType &shape)
SizeType m_OffSetTable[Dimension]
Definition: utlNDArray.h:1167
NDArray(const NDArray< TMatrixValueType, 4 > &r)
Superclass::Pointer Pointer
Base class for utl::NDArray.
Definition: utlMatrix.h:22
Superclass::ScalarValueType ScalarValueType
bool IsMinorSymmetric() const
Definition: utl.h:90
Superclass::SizeType SizeType
Definition: utlNDArray.h:154
void Fill(const T &value)
Definition: utlNDArray.h:922
SizeType GetOffset(const ShapeType &shapeIndex) const
Definition: utlNDArray.h:338
NDArray(const NDArray< T, 4 > &mat)
NDArray< T, 4 > & operator=(NDArray< T, 4 > &&r)
void Swap(NDArrayBase< T, Dim > &vec)
Definition: utlNDArray.h:1062
NDArray(NDArray< T, 4 > &&r)
Superclass::SizeType SizeType
bool IsMajorSymmetric() const
const ValueType * ConstPointer
Definition: utlNDArray.h:175
NDArray< T, 2 > GetRefSubMatrix(const int i, const int j) const
NDArray(const ShapeType &shape, const T r)
bool IsSame(const T &value, const T &v0, const double eps=1e-10)
Definition: utlCore.h:1337
void EigenDecompositionNonSymmetricMatrix(NDArray< T, 1 > &valReal, NDArray< T, 1 > &valImg) const
Definition: utlMatrix.h:408
Superclass::ShapeType ShapeType
void Convert1To2Tensor(const utl::NDArray< T, 1 > &vec, utl::NDArray< T, 2 > &mat)
Definition: utlDMRI.h:386
void EigenDecompositionWithMinorSymmetry(NDArray< T, 1 > &valReal, NDArray< T, 1 > &valImg, std::vector< NDArray< T, 2 > > &matRealR, std::vector< NDArray< T, 2 > > &matImgR, std::vector< NDArray< T, 2 > > &matRealL, std::vector< NDArray< T, 2 > > &matImgL) const
SizeType Size() const
Definition: utlNDArray.h:321
Superclass::Iterator Iterator
NDArrayBase< T, Dim > & operator=(const Expr< EType, typename EType::ValueType > &src)
Definition: utlNDArray.h:385
NDArray(const T *vec, const ShapeType &shape)
utl::remove_complex_t< T > ScalarValueType
Definition: utlNDArray.h:165
void SetData(T *const data, const unsigned rows, const unsigned cols)
Definition: utlMatrix.h:190
Superclass::ConstIterator ConstIterator
Superclass::ValueType ValueType
Superclass::ShapeType ShapeType
Definition: utlNDArray.h:155
NDArray< T, 4 > & operator=(NDArray< T, 4 > &r)
#define __utl_ndarray_alloc_blah(shape)
Definition: utlNDArray.h:30
void CopyData(T *const data, const ShapeType &shape)
Definition: utlNDArray.h:689
NDArray<T,4> is a 4th order tensor.
NDArrayBase< T, 4 > Superclass