DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
utlDMRI.h
Go to the documentation of this file.
1 
10 #ifndef __utlDMRI_h
11 #define __utlDMRI_h
12 
13 #include "utlCore.h"
14 #include "utlNDArray.h"
15 #include "utlSTDHeaders.h"
16 
20 enum {
25 };
26 
27 enum {
32 };
33 
34 enum {
37 };
38 
39 enum {
42 };
43 
44 // enum {
45 // CARTESIAN=0,
46 // SPHERICAL=1
47 // };
48 
49 enum{
50  // 9d [xx, xy, xz, yx, yy, yz, zx, zy, zz]
52  // [xx, xy, xz, yy, yz, zz], fsl
54  // [xx, yx, yy, zx, zy, zz]
56  // [xx, yy, zz, xy, xz yz], mrtrix, internal in vistasoft, AFQ
58  // [xx, yy, zz, sqrt(2)*xy, sqrt(2)*xz, sqrt(2)*yz], embed 3x3 matrix into 6x1 vector
60 };
61 
62 enum{
67 };
68 
69 enum{
72 };
73 
74 namespace utl
75 {
76 
77 
78 // [>* The name and dimention map for vectors. <]
79 // typedef utl_unordered_map<std::string, int> VecDimType;
80 // static VecDimType vectorNameDimMap{
81 // {"scalar",1},
82 // {"frame", 9},
83 // {"point",3}
84 // };
85 
92 inline int
93 GetScalarsDimentionByName(const std::string& name)
94 {
95  // auto ii = vectorNameDimMap.find(name);
96  // if (ii!=vectorNameDimMap.end())
97  // return ii->second;
98  // else
99  // return 1;
100 
101  size_t found = name.find_last_of("_");
102  if (found==std::string::npos)
103  return 1;
104 
105  std::string dim = name.substr(found+1);
106  if (utl::IsInt(dim))
107  return utl::ConvertStringToNumber<int>(dim);
108  else
109  return 1;
110 }
111 
113 template <class T>
114 inline std::vector<T>
115 GetScalarsByName(const std::vector<T>& vec, const std::vector<std::string>& nameVec, const std::string& name)
116 {
117  int start=0;
118  std::vector<T> result;
119  for ( int i = 0; i < nameVec.size(); ++i )
120  {
121  int len = GetScalarsDimentionByName(nameVec[i]);
122  if (nameVec[i]==name)
123  {
124  for ( int j = 0; j < len; ++j )
125  result.push_back(vec[start+j]);
126  return result;
127  }
128  start += len;
129  }
130  return result;
131 }
132 
133 template <class T>
134 inline void
135 SetScalarsByName(const std::vector<T>& vec, const std::vector<std::string>& nameVec, const std::vector<T>& scalars, const std::string& name)
136 {
137  int start=0;
138  for ( int i = 0; i < nameVec.size(); ++i )
139  {
140  int len = GetScalarsDimentionByName(nameVec[i]);
141  utlSAException(len!=scalars.size())(nameVec[i])(len)(scalars.size()).msg("the size of scalars is not consistent");
142  if (nameVec[i]==name)
143  {
144  for ( int j = 0; j < len; ++j )
145  vec[start+j] = scalars[j];
146  }
147  start += len;
148  }
149 }
150 
151 template <class T>
152 inline void
153 RemoveScalarsByName(std::vector<T>& vec, const std::vector<std::string>& nameVec, const std::string& name)
154 {
155  int start=0, len=0, i;
156  for ( i = 0; i < nameVec.size(); ++i )
157  {
158  len = GetScalarsDimentionByName(nameVec[i]);
159  if (nameVec[i]==name)
160  break;
161  start += len;
162  }
163 
164  // no name found
165  if (i==nameVec.size())
166  return;
167 
168  // // remove name
169  // int jj = i;
170  // nameVec.erase(nameVec.begin()+jj);
171 
172  // remove scalars
173  for ( i = start; i+len < vec.size(); ++i )
174  vec[i] = vec[i+len];
175 
176  vec.resize(vec.size()-len);
177 }
178 
179 
181 inline int
182 DimToRankSH ( const int dimm )
183 {
184  double result = -1.5 + std::sqrt(0.25 + 2.0*dimm);
185 
186  utlGlobalAssert(IsInt(result) && IsEven(result), "Logical ERROR! Rank of SH here must be an even integer. Wrong dimension! rank=" << result <<", dimm = " << dimm);
187 
188  return (int) result;
189 }
190 
192 inline int
193 RankToDimSH (const int shRank)
194 {
195  utlSAGlobalException(!utl::IsEven(shRank))(shRank).msg("shRank should be even");
196  if (shRank<0)
197  return 0;
198  else
199  return (shRank + 1)*(shRank + 2)/2;
200 }
201 
203 inline int
204 GetIndexSHj(const int l, const int m)
205 {
206  return (l*l + l)/2 + m;
207 }
208 
210 inline std::vector<int>
211 GetIndexSHlm(const int j)
212 {
213  std::vector<int> lm(2,0);
214 
215  if (j==0)
216  return lm;
217 
218  int l = std::ceil(0.5*(std::sqrt(9.+8.*j)-3.0));
219  if (utl::IsOdd(l))
220  l+=1;
221  int residual = j+1 - (l-1)*l/2;
222  lm[0] = l;
223  lm[1] = residual-1 - l;
224 
225  return lm;
226 }
227 
233 template < typename T >
234 std::vector<T>
235 GetE1E2FromFAMD ( const T fa, const T meanEigenValue, const bool isE2E3Equal=true )
236 {
237  std::vector<T> eigenValues(2);
238  double tmp = fa/std::sqrt(3.0-2.0*fa*fa);
239  if (isE2E3Equal)
240  {
241  eigenValues[0] = meanEigenValue*(1.0 + 2.0*tmp);
242  eigenValues[1] = meanEigenValue*(1.0 - tmp);
243  }
244  else
245  {
246  utlSAException(fa>1.0/utl::SQRT2)(fa).msg("fa is too large so that e3 is negative");
247  eigenValues[0] = meanEigenValue*(1.0 + tmp);
248  eigenValues[1] = meanEigenValue*(1.0 - 2.0*tmp);
249  }
250  return eigenValues;
251 }
252 
255 template<typename V1Type, typename V2Type>
256 void
257 ConvertTensor6DTo9D(const V1Type& v6d, V2Type& v9d, int v6dStoreWay)
258 {
259  if (v6dStoreWay==TENSOR_LOWER_TRIANGULAR)
260  {
261  v9d[0] = v6d[0]; // xx
262  v9d[1] = v6d[1]; // xy
263  v9d[2] = v6d[3]; // xz
264  v9d[3] = v6d[1]; // yx
265  v9d[4] = v6d[2]; // yy
266  v9d[5] = v6d[4]; // yz
267  v9d[6] = v6d[3]; // zx
268  v9d[7] = v6d[4]; // zy
269  v9d[8] = v6d[5]; // zz
270  }
271  else if (v6dStoreWay==TENSOR_UPPER_TRIANGULAR)
272  {
273  v9d[0] = v6d[0]; // xx
274  v9d[1] = v6d[1]; // xy
275  v9d[2] = v6d[2]; // xz
276  v9d[3] = v6d[1]; // yx
277  v9d[4] = v6d[3]; // yy
278  v9d[5] = v6d[4]; // yz
279  v9d[6] = v6d[2]; // zx
280  v9d[7] = v6d[4]; // zy
281  v9d[8] = v6d[5]; // zz
282  }
283  else if (v6dStoreWay==TENSOR_EMBED6D)
284  {
285  v9d[0] = v6d[0]; // xx
286  v9d[1] = v6d[3]*utl::SQRT1_2; // xy
287  v9d[2] = v6d[4]*utl::SQRT1_2; // xz
288  v9d[3] = v6d[3]*utl::SQRT1_2; // yx
289  v9d[4] = v6d[1]; // yy
290  v9d[5] = v6d[5]*utl::SQRT1_2; // yz
291  v9d[6] = v6d[4]*utl::SQRT1_2; // zx
292  v9d[7] = v6d[5]*utl::SQRT1_2; // zy
293  v9d[8] = v6d[2]; // zz
294  }
295  else if (v6dStoreWay==TENSOR_DIAGONAL_FIRST)
296  {
297  v9d[0] = v6d[0]; // xx
298  v9d[1] = v6d[3]; // xy
299  v9d[2] = v6d[4]; // xz
300  v9d[3] = v6d[3]; // yx
301  v9d[4] = v6d[1]; // yy
302  v9d[5] = v6d[5]; // yz
303  v9d[6] = v6d[4]; // zx
304  v9d[7] = v6d[5]; // zy
305  v9d[8] = v6d[2]; // zz
306  }
307  else
308  utlException(true, "wrong v1StoreWay");
309 }
310 
313 template<typename V1Type, typename V2Type>
314 void
315 ConvertTensor9DTo6D(const V1Type& v9d, V2Type& v6d, int v6dStoreWay)
316 {
317  if (v6dStoreWay==TENSOR_LOWER_TRIANGULAR)
318  {
319  v6d[0] = v9d[0]; // xx
320  v6d[1] = v9d[3]; // yx
321  v6d[2] = v9d[4]; // yy
322  v6d[3] = v9d[6]; // zx
323  v6d[4] = v9d[7]; // zy
324  v6d[5] = v9d[8]; // zz
325  }
326  else if (v6dStoreWay==TENSOR_UPPER_TRIANGULAR)
327  {
328  v6d[0] = v9d[0]; // xx
329  v6d[1] = v9d[1]; // xy
330  v6d[2] = v9d[2]; // xz
331  v6d[3] = v9d[4]; // yy
332  v6d[4] = v9d[5]; // yz
333  v6d[5] = v9d[8]; // zz
334  }
335  else if (v6dStoreWay==TENSOR_EMBED6D)
336  {
337  v6d[0] = v9d[0]; // xx
338  v6d[1] = v9d[4]; // yy
339  v6d[2] = v9d[8]; // zz
340  v6d[3] = v9d[1]*utl::SQRT2; // sqrt(2)*xy
341  v6d[4] = v9d[2]*utl::SQRT2; // sqrt(2)*xz
342  v6d[5] = v9d[5]*utl::SQRT2; // sqrt(2)*yz
343  }
344  else if (v6dStoreWay==TENSOR_DIAGONAL_FIRST)
345  {
346  v6d[0] = v9d[0]; // xx
347  v6d[1] = v9d[4]; // yy
348  v6d[2] = v9d[8]; // zz
349  v6d[3] = v9d[1]; // xy
350  v6d[4] = v9d[2]; // xz
351  v6d[5] = v9d[5]; // yz
352  }
353  else
354  utlException(true, "wrong v1StoreWay");
355 }
356 
357 template<typename V1Type, typename V2Type>
358 void
359 ConvertTensor6DTo6D(const V1Type& v6d1, V2Type& v6d2, int s1, int s2)
360 {
361  if (s1==s2)
362  {
363  for ( int i = 0; i < 6; ++i )
364  v6d2[i] = v6d1[i];
365  }
366  else
367  {
368  std::vector<double> v9d(9);
369  ConvertTensor6DTo9D(v6d1, v9d, s1);
370  ConvertTensor9DTo6D(v9d, v6d2, s2);
371  }
372 }
373 
375 template <class T>
376 inline void
378 {
379  utlException(mat.Rows()!=3 || mat.Cols()!=3, "It should be in 3D space");
380  vec.ReSize(6);
382 }
383 
384 template <class T>
385 inline void
387 {
388  utlException(vec.Size()!=6, "It should be in 3D space");
389  mat.ReSize(3,3);
391 }
392 
393 template <class T>
394 inline void
396 {
397  if (grad.Size()==0)
398  return;
399  utlSAException(grad.Cols()!=3)(grad.Rows())(grad.Cols())("grad should be Nx3 matrix");
400  utl::NDArray<T,1> vec;
401  for ( int i = 0; i < grad.Rows(); ++i )
402  {
403  vec = grad.GetRow(i);
404  double norm = vec.GetTwoNorm();
405  if (norm>1e-10)
406  {
407  vec /= norm;
408  grad.SetRow(i, vec);
409  }
410  }
411 }
412 
413 template <class T>
414 inline utl_shared_ptr<NDArray<T,2> >
415 ReadGrad(const std::string& grad_str, const int NoSymmetricDuple=DIRECTION_NODUPLICATE, const int mode=CARTESIAN_TO_SPHERICAL,
416  const int flipx=DIRECTION_NOFLIP, const int flipy=DIRECTION_NOFLIP, const int flipz=DIRECTION_NOFLIP, const bool need_normalize=true)
417 {
418  std::vector < std::vector<std::string> > strVec;
419 
420  // use three delimiters ' ', ',', '\t'
421  ReadLinesFirstlineCheck(grad_str, strVec, " \t,");
422 
423  for ( int i = 0; i < strVec.size(); i += 1 )
424  utlGlobalException(strVec[i].size()!=3, "wrong dimension in gradient file! strVec[i].size()="<< strVec[i].size() <<", grad_str="<<grad_str);
425 
426  utl_shared_ptr<NDArray<T,2> > grad (new NDArray<T,2>() );
427  if (NoSymmetricDuple==DIRECTION_DUPLICATE)
428  grad->ReSize(2*strVec.size(),3);
429  else
430  grad->ReSize(strVec.size(),3);
431 
432  std::vector<double> xyz(3);
433  int jj=0;
434  for (int i=0; i<strVec.size(); i++)
435  {
436  utlGlobalException(strVec[i].size()!=3, "wrong size for grad");
437  std::istringstream ( strVec[i][0] ) >> xyz[0];
438  std::istringstream ( strVec[i][1] ) >> xyz[1];
439  std::istringstream ( strVec[i][2] ) >> xyz[2];
440 
442  {
443  if (need_normalize)
444  xyz = NormalizeUnitNorm(xyz);
445  xyz[0] = (flipx==DIRECTION_FLIP?-1:1)*xyz[0]; xyz[1] = (flipy==DIRECTION_FLIP?-1:1)*xyz[1]; xyz[2] = (flipz==DIRECTION_FLIP?-1:1)*xyz[2];
446  if (mode==CARTESIAN_TO_SPHERICAL)
447  cartesian2Spherical(xyz[0],xyz[1],xyz[2]);
448  }
449  else if(mode==SPHERICAL_TO_CARTESIAN || mode==SPHERICAL_TO_SPHERICAL)
450  {
451  spherical2Cartesian(xyz[0],xyz[1],xyz[2]);
452  if (need_normalize)
453  xyz = NormalizeUnitNorm(xyz);
454  xyz[0] = (flipx==DIRECTION_FLIP?-1:1)*xyz[0]; xyz[1] = (flipy==DIRECTION_FLIP?-1:1)*xyz[1]; xyz[2] = (flipz==DIRECTION_FLIP?-1:1)*xyz[2];
455  if (mode==SPHERICAL_TO_SPHERICAL)
456  cartesian2Spherical(xyz[0],xyz[1],xyz[2]);
457  }
458  else
459  utlGlobalException(true, "wrong mode");
460 
461  (*grad)(jj,0)=xyz[0], (*grad)(jj,1)=xyz[1], (*grad)(jj,2)=xyz[2];
462  jj++;
463 
464  if(NoSymmetricDuple==DIRECTION_DUPLICATE)
465  {
467  { xyz[0] = -xyz[0]; xyz[1] = -xyz[1]; xyz[2] = -xyz[2]; }
468  else
469  { xyz[1]=M_PI-xyz[1]; xyz[2]=M_PI+xyz[2]; xyz[2]=(xyz[2]>=2*M_PI)?(xyz[2]-2*M_PI):xyz[2]; }
470  (*grad)(jj,0)=xyz[0], (*grad)(jj,1)=xyz[1], (*grad)(jj,2)=xyz[2];
471  jj++;
472  }
473  }
474 
475  return grad;
476 }
477 
478 template < class T >
479 void
480 MatchBVectorAndGradientMatrix (const T& br, std::vector<T>& vec, const NDArray<T,2>& grad )
481 {
482  utlException(grad.Rows()==0, "grad.size()==0");
483  vec = std::vector<T>(grad.Rows(), br);
484 }
485 
486 template < class T >
487 void
488 MatchBVectorAndGradientMatrix ( std::vector<T>& vec, NDArray<T,2>& grad )
489 {
490  utlException(vec.size()==0 || grad.Rows()==0, "wrong size! vec.size()="<<vec.size()<<", grad.Rows()="<<grad.Rows());
491 
492  if (grad.Rows()==vec.size())
493  return;
494 
495  std::vector<T> vec_result;
496  if (vec.size()==1)
497  {
498  MatchBVectorAndGradientMatrix(vec[0], vec_result, grad);
499  vec = vec_result;
500  return;
501  }
502 
503  vec_result.resize(vec.size()* grad.Rows());
504  NDArray<T,2> grad_result(vec.size()* grad.Rows(), grad.Columns());
505  for (int i = 0; i<vec.size(); i++)
506  {
507  for ( int j = 0; j < grad.Rows(); j += 1 )
508  {
509  vec_result[i*grad.Rows()+j] = vec[i];
510  for ( int kk = 0; kk < 3; kk += 1 )
511  grad_result(i*grad.Rows()+j,kk) = grad(j,kk);
512  }
513  }
514  vec = vec_result;
515  grad = grad_result;
516  return;
517 }
518 
519 inline int
520 GetFiberTractsFormatFromFileExtension(const std::string& filename)
521 {
522  std::string ext, fileNoExt;
523  utl::GetFileExtension(filename, ext, fileNoExt);
524  int format = TRACTS_UNKNOWN;
525  if (ext=="trk")
526  format = TRACTS_TRK;
527  else if (ext=="tck")
528  format = TRACTS_TCK;
529  else if (ext=="vtk")
530  format = TRACTS_VTK;
531  else
532  format = TRACTS_UNKNOWN;
533  return format;
534 }
535 
536 inline int
537 GetFiberTractsFormat(const std::string& filename)
538 {
539  std::string ext, fileNoExt;
540  utl::GetFileExtension(filename, ext, fileNoExt);
541 
542  FILE* file;
543  file = fopen(filename.c_str(), "rb");
544  if (!file)
545  utlGlobalException(true, "Unable to open file " + filename);
546 
547  int format = TRACTS_UNKNOWN;
548  if (ext=="trk")
549  {
550  format = TRACTS_TRK;
551  char code[6];
552  fseek (file, 0, SEEK_SET);
553  size_t rsize = fread(code, 1, 6, file);
554  utlGlobalException(strcmp(code,"TRACK")!=0, "Wrong data format TRACTS_TRK. Magic code is wrong. It is " + std::string(code) +". Should be 'TRACK'");
555  utlGlobalException(rsize!=6, "wrong size read");
556  }
557  else if (ext=="tck")
558  {
559  format = TRACTS_TCK;
560  }
561  else if (ext=="vtk")
562  {
563  format = TRACTS_VTK;
564  }
565  else
566  format = TRACTS_UNKNOWN;
567  fclose(file);
568  return format;
569 }
570 
581 template <class T>
583 GetDTIDesignMatrix ( const utl::Matrix<T>& gradMat, const std::vector<T>& bVec, int dwi_normalize )
584 {
585  utlSAException(gradMat.Rows()!=bVec.size())(gradMat.Rows())(bVec.size()).msg("wrong size of gradMat and bVec");
586  utl::Matrix<T> mat;
587  if (dwi_normalize==DWI_NORMALIZE)
588  mat.ReSize(gradMat.Rows(), 6);
589  else if (dwi_normalize==DWI_NONORMALIZE)
590  mat.ReSize(gradMat.Rows(), 7);
591  else
592  utlException(true, "wrong dwi_normalize");
593 
594  for ( int i = 0; i < gradMat.Rows(); ++i )
595  {
596  mat(i,0) = -1.0 * bVec[i] * gradMat(i,0) * gradMat(i,0); // xx
597  mat(i,1) = -2.0 * bVec[i] * gradMat(i,0) * gradMat(i,1); // xy
598  mat(i,2) = -2.0 * bVec[i] * gradMat(i,0) * gradMat(i,2); // xz
599  mat(i,3) = -1.0 * bVec[i] * gradMat(i,1) * gradMat(i,1); // yy
600  mat(i,4) = -2.0 * bVec[i] * gradMat(i,1) * gradMat(i,2); // yz
601  mat(i,5) = -1.0 * bVec[i] * gradMat(i,2) * gradMat(i,2); // zz
602  if (dwi_normalize==DWI_NONORMALIZE)
603  mat(i,6) = 1; // S0
604  }
605  return mat;
606 }
607 
608 /* Function taken from 3D Slicer, vtkDiffusionTensorMathematics
609  *
610  * This is sort of the inverse of code from Gordon Kindlmann for mapping
611  * the mode (index value) to RGB. See vtkTensorMathematics for that code.
612  * There may be a simpler way to do this but this works.
613  * Note this expects a "0 1" Hue Range in the vtkLookupTable used to
614  * display the glyphs.
615  */
616 inline void
617 RGBToIndex(double R, double G, double B, double &index)
618 {
619 
620  // remove the gray part of the color.
621  // this is so we can use the model where either R,G, or B is 0.
622  // then we scale so that the max of the other two is one.
623  double min = R;
624  int minIdx = 0;
625 
626  if (G < min)
627  {
628  min = G;
629  minIdx = 1;
630  }
631  if (B < min)
632  {
633  min = B;
634  minIdx = 2;
635  }
636 
637  // make the smallest of R,G,B equal 0
638  R = R - min;
639  G = G - min;
640  B = B - min;
641 
642  // now take the max, and scale it to be 1.
643  double max = R;
644  int maxIdx = 0;
645  if (G > max)
646  {
647  max = G;
648  maxIdx = 1;
649  }
650  if (B > max)
651  {
652  max = B;
653  maxIdx = 2;
654  }
655 
656  R = R/max;
657  G = G/max;
658  B = B/max;
659 
660 
661  // now using the inverse sextants, map this into an index.
662  // switch (sextant) {
663  // case 0: { R = 1; G = frac; B = 0; break; }
664  // case 1: { R = 1-frac; G = 1; B = 0; break; }
665  // case 2: { R = 0; G = 1; B = frac; break; }
666  // case 3: { R = 0; G = 1-frac; B = 1; break; }
667  // case 4: { R = frac; G = 0; B = 1; break; }
668  // case 5: { R = 1; G = 0; B = 1-frac; break; }
669  // }
670  int sextant;
671  if (maxIdx == 0 && minIdx == 2) sextant = 0;
672  if (maxIdx == 1 && minIdx == 2) sextant = 1;
673  if (maxIdx == 1 && minIdx == 0) sextant = 2;
674  if (maxIdx == 2 && minIdx == 0) sextant = 3;
675  if (maxIdx == 2 && minIdx == 1) sextant = 4;
676  if (maxIdx == 0 && minIdx == 1) sextant = 5;
677 
678  double offset;
679  offset = 256/6;
680 
681  switch (sextant)
682  {
683  case 0: { index = G*offset; break; }
684  case 1: { index = offset + (1-R)*offset; break; }
685  case 2: { index = offset*2 + B*offset; break; }
686  case 3: { index = offset*3 + (1-G)*offset; break; }
687  case 4: { index = offset*4 + R*offset; break; }
688  case 5: { index = offset*5 + (1-B)*offset; break; }
689  }
690 }
691 
692 }
693 
696 #endif
int GetFiberTractsFormatFromFileExtension(const std::string &filename)
Definition: utlDMRI.h:520
void NormalizeGrad(const utl::NDArray< T, 2 > &grad)
Definition: utlDMRI.h:395
NDArray<T,1> is a vector class which uses blas mkl.
Definition: utlVector.h:36
std::vector< T > GetE1E2FromFAMD(const T fa, const T meanEigenValue, const bool isE2E3Equal=true)
Definition: utlDMRI.h:235
bool IsOdd(const int value)
Definition: utlCore.h:825
std::shared_ptr< NDArray< T, 2 > > ReadGrad(const std::string &grad_str, 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)
Definition: utlDMRI.h:415
void GetRow(const int index, T *vec) const
Definition: utlMatrix.h:266
int GetIndexSHj(const int l, const int m)
Definition: utlDMRI.h:204
bool IsEven(const int value)
Definition: utlCore.h:819
void MatchBVectorAndGradientMatrix(const T &br, std::vector< T > &vec, const NDArray< T, 2 > &grad)
Definition: utlDMRI.h:480
static constexpr double SQRT1_2
Definition: utlConstants.h:33
SizeType Columns() const
Definition: utlMatrix.h:137
void SetScalarsByName(const std::vector< T > &vec, const std::vector< std::string > &nameVec, const std::vector< T > &scalars, const std::string &name)
Definition: utlDMRI.h:135
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
void spherical2Cartesian(const T r, const T theta, const T phi, T &x, T &y, T &z)
Definition: utlCore.h:1602
std::vector< int > GetIndexSHlm(const int j)
Definition: utlDMRI.h:211
SizeType Rows() const
Definition: utlMatrix.h:136
#define utlSAGlobalException(expr)
Definition: utlCoreMacro.h:362
NDArray<T,2> is a row-major matrix.
Definition: utlMatrix.h:37
void ReadLinesFirstlineCheck(const std::string &filename, std::vector< std::vector< std::string > > &strVec, const char *cc=" ")
Definition: utlCore.h:890
void RemoveScalarsByName(std::vector< T > &vec, const std::vector< std::string > &nameVec, const std::string &name)
Definition: utlDMRI.h:153
utl::Matrix< T > GetDTIDesignMatrix(const utl::Matrix< T > &gradMat, const std::vector< T > &bVec, int dwi_normalize)
Definition: utlDMRI.h:583
int DimToRankSH(const int dimm)
Definition: utlDMRI.h:182
const T & max(const T &a, const T &b)
Return the maximum between a and b.
Definition: utlCore.h:263
double GetTwoNorm() const
Definition: utlNDArray.h:1013
NDArray< T, 2 > & SetRow(const int index, T const *vec)
Definition: utlMatrix.h:241
int GetFiberTractsFormat(const std::string &filename)
Definition: utlDMRI.h:537
const T & min(const T &a, const T &b)
Return the minimum between a and b.
Definition: utlCore.h:257
void Convert2To1Tensor(const utl::NDArray< T, 2 > &mat, utl::NDArray< T, 1 > &vec)
Definition: utlDMRI.h:377
#define M_PI
Definition: utlCoreMacro.h:57
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
#define utlGlobalAssert(cond, expout)
Definition: utlCoreMacro.h:381
Definition: utl.h:90
bool IsInt(const std::string &input, const double epss=1e-10)
Definition: utlCore.h:792
SizeType Cols() const
Definition: utlMatrix.h:138
#define utlSAException(expr)
Definition: utlCoreMacro.h:543
int RankToDimSH(const int shRank)
Definition: utlDMRI.h:193
void ConvertTensor6DTo9D(const V1Type &v6d, V2Type &v9d, int v6dStoreWay)
Definition: utlDMRI.h:257
std::vector< T > GetScalarsByName(const std::vector< T > &vec, const std::vector< std::string > &nameVec, const std::string &name)
Definition: utlDMRI.h:115
void RGBToIndex(double R, double G, double B, double &index)
Definition: utlDMRI.h:617
int GetScalarsDimentionByName(const std::string &name)
Definition: utlDMRI.h:93
void ConvertTensor6DTo6D(const V1Type &v6d1, V2Type &v6d2, int s1, int s2)
Definition: utlDMRI.h:359
void cartesian2Spherical(const T x, const T y, const T z, T &r, T &theta, T &phi)
Definition: utlCore.h:1579
VectorType NormalizeUnitNorm(const VectorType &v, const int nSize)
Definition: utlCore.h:672
void Convert1To2Tensor(const utl::NDArray< T, 1 > &vec, utl::NDArray< T, 2 > &mat)
Definition: utlDMRI.h:386
SizeType Size() const
Definition: utlNDArray.h:321
bool ReSize(const SizeType rows, const SizeType cols)
Definition: utlMatrix.h:140
bool ReSize(const SizeType size)
Definition: utlVector.h:167
void GetFileExtension(const std::string &fileNameAbsolute, std::string &ext, std::string &fileNoExt)
Definition: utlCore.h:559
void ConvertTensor9DTo6D(const V1Type &v9d, V2Type &v6d, int v6dStoreWay)
Definition: utlDMRI.h:315
static constexpr double SQRT2
Definition: utlConstants.h:31