DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkSamplingSchemeQSpaceIncrementalEstimationFilter.hxx
Go to the documentation of this file.
1 
18 #ifndef __itkSamplingSchemeQSpaceIncrementalEstimationFilter_hxx
19 #define __itkSamplingSchemeQSpaceIncrementalEstimationFilter_hxx
20 
21 #include "itkSamplingSchemeQSpaceIncrementalEstimationFilter.h"
22 #include "utl.h"
23 
24 namespace itk
25 {
26 
27 template< class TSamplingType >
28 SamplingSchemeQSpaceIncrementalEstimationFilter< TSamplingType >
29 ::SamplingSchemeQSpaceIncrementalEstimationFilter()
30  : Superclass(),
31  m_FineOrientations(new MatrixType())
32 {
33  m_TessellationOrder = 7;
34 }
35 
36 template< class TSamplingType >
37 void
38 SamplingSchemeQSpaceIncrementalEstimationFilter< TSamplingType >
39 ::Initialization()
40 {
41  // SamplingType* input = const_cast< SamplingType * >( this->GetInput() );
42  if (this->m_FineOrientations->Size()==0)
43  {
44  this->m_FineOrientations = utl::ReadGrad<double>(this->m_TessellationOrder, DIRECTION_NODUPLICATE, CARTESIAN_TO_CARTESIAN); // catesian
45  }
46 
47  if (this->m_InitialOrientations)
48  {
49  MatrixPointer initialMatrix = this->m_InitialOrientations->GetOrientationsCartesian();
50  *this->m_FineOrientations = utl::ConnectUtlMatrix(*initialMatrix, *this->m_FineOrientations, true);
51 
52  SamplingPointer output = this->GetOutputOrientations();
53  output->SetOrientationsCartesian(this->m_FineOrientations);
54  output->SetIndicesInShells(this->m_InitialOrientations->GetIndicesInShells());
55  Index2DVectorPointer indices = output->GetIndicesInShells();
56 
57  unsigned int numberInitialSamples = this->m_InitialOrientations->GetNumberOfSamples();
58  unsigned int num=0;
59  for ( unsigned int i = 0; i < indices->size(); i += 1 )
60  {
61  num += indices->operator[](i).size();
62  }
63  utlGlobalException(num!=numberInitialSamples, "wrong intialization in m_InitialOrientations");
64  }
65 }
66 
67 // template< class TSamplingType >
68 // void
69 // SamplingSchemeQSpaceIncrementalEstimationFilter< TSamplingType >
70 // ::IndicesOfMaximalDistance(const MatrixType& mat, unsigned int& r1, unsigned int& r2)
71 // {
72 // unsigned int row = mat.Rows();
73 // double x,y,z,x1,y1,z1, dotMin=100, dot;
74 // for ( unsigned int i = 1; i < row; i += 1 )
75 // {
76 // x=mat(i,0), y=mat(i,1), z=mat(i,2);
77 // for ( unsigned int j = 0; j < i; j += 1 )
78 // {
79 // x1=mat(i,0), y1=mat(i,1), z1=mat(i,2);
80 // dot = x*x1 + y*y1 + z*z1;
81 // if ( dot<0 )
82 // dot = -dot;
83 
84 // if (dot<dotMin)
85 // {
86 // dotMin = dot;
87 // r1=i, r2=j;
88 // }
89 // }
90 // }
91 // }
92 
93 template< class TSamplingType >
94 void
95 SamplingSchemeQSpaceIncrementalEstimationFilter< TSamplingType >
96 ::GenerateData()
97 {
98  Initialization();
99 
100  unsigned int numberOfShells = this->m_NumbersInShell.size();
101  unsigned int numberInitialSamples = 0;
102  if (this->m_InitialOrientations)
103  numberInitialSamples = this->m_InitialOrientations->GetNumberOfSamples();
104 
105  SamplingPointer output = this->GetOutputOrientations();
106  output->SetOrientationsCartesian(this->m_FineOrientations);
107  Index2DVectorPointer indices = output->GetIndicesInShells();
108 
109  unsigned int numberOfSamples = output->GetNumberOfSamples();
110  MatrixPointer outputMatrix = output->GetOrientationsCartesian();
111 
112  // itk::TimeProbe clock;
113  // clock.Start();
114  // MatrixPointer dotMatrix (new MatrixType());
115  // MatrixPointer elecMatrix (new MatrixType());
116  // if (numberOfShells>1)
117  // {
118  // if (this->m_CriteriaType==Self::DISTANCE)
119  // dotMatrix = output->CalculateInnerProductMatrix(true);
120  // else if (this->m_CriteriaType==Self::ELECTROSTATIC)
121  // elecMatrix = output->CalculateElectrostaticEnergyMatrix(this->m_ElectrostaticOrder);
122  // }
123  // clock.Stop();
124  // std::cout << "after matrix calculation: " << clock.GetTotal() << std::endl;
125 
126  utlGlobalException(this->m_NumbersInShell.size()==0, "need to set this->m_NumbersInShell");
127  unsigned int totalNumbers=0;
128  for ( unsigned int i = 0; i < this->m_NumbersInShell.size(); i += 1 )
129  {
130  totalNumbers += this->m_NumbersInShell[i];
131  }
132 
133  std::vector<int> chosenIndex;
134  std::vector<char> hasChosen(numberOfSamples, 0);
135  std::cout << "numberOfShells = " << numberOfShells << std::endl << std::flush;
136  std::cout << "numberOfSamples = " << numberOfSamples << std::endl << std::flush;
137 
138  if (!this->m_InitialOrientations)
139  {
140  // first two points
141  int d0, d1;
142  // clock.Start();
143  // if (numberOfShells>1)
144  // {
145  // if (this->m_CriteriaType==Self::DISTANCE)
146  // utl::ArgminSymmetricMatrix<MatrixType>(*dotMatrix, numberOfSamples, d0, d1, false);
147  // else if (this->m_CriteriaType==Self::ELECTROSTATIC)
148  // utl::ArgminSymmetricMatrix<MatrixType>(*elecMatrix, numberOfSamples, d0, d1, false);
149  // }
150 
151  IndexVectorType indexTmp;
152  // utlPrintVar2(true, d0, d1);
153  // clock.Stop();
154  // std::cout << "after matrix minimal: " << clock.GetTotal() << std::endl;
155  if (numberOfShells==1)
156  {
157  d0=0;
158  indexTmp.push_back(d0);
159  // indexTmp.push_back(d1);
160  indices->push_back(indexTmp);
161  }
162  else
163  {
164  d0=0;
165  indexTmp.push_back(d0);
166  indices->push_back(indexTmp);
167  // indexTmp[0] = d1;
168  // indices->push_back(indexTmp);
169  }
170  chosenIndex.push_back(d0);
171  // chosenIndex.push_back(d1);
172  hasChosen[d0] = 1;
173  // hasChosen[d1] = 1;
174  }
175 
176  for ( unsigned int i = 0; i < indices->size(); i += 1 )
177  for ( unsigned int j = 0; j < (*indices)[i].size(); j += 1 )
178  {
179  int ii = (*indices)[i][j];
180  chosenIndex.push_back(ii);
181  hasChosen[ii] = 1;
182  }
183 
184  if (numberOfShells==1)
185  {
186 
187  if (this->m_CriteriaType==Self::DISTANCE)
188  {
189  double x,y,z,x1,y1,z1;
190  std::vector<double> maxDots(numberOfSamples, -1);
191  // clock.Start();
192  for ( unsigned int i = (this->m_InitialOrientations?numberInitialSamples:1); i < totalNumbers; i += 1 )
193  {
194  int index=-1;
195  double minValue=std::numeric_limits<double>::max();
196  for ( unsigned int j = 0; j < numberOfSamples; j += 1 )
197  {
198  if (hasChosen[j]==0)
199  {
200  double maxValue_j = -std::numeric_limits<double>::max();
201  x = output->operator[](j)[0];
202  y = output->operator[](j)[1];
203  z = output->operator[](j)[2];
204  if (maxDots[j]< 0 )
205  {
206  for ( unsigned int k = 0; k < chosenIndex.size(); k += 1 )
207  {
208  unsigned int kk = chosenIndex[k];
209  // double value = dotMatrix->operator()(j,kk);
210  x1 = output->operator[](kk)[0];
211  y1 = output->operator[](kk)[1];
212  z1 = output->operator[](kk)[2];
213  double value = x*x1 + y*y1 + z*z1;
214  if (value<0)
215  value = -value;
216  if (value>maxValue_j)
217  maxValue_j = value;
218  }
219  maxDots[j] = maxValue_j;
220  }
221  else
222  {
223  unsigned int kk = chosenIndex.back();
224  // double value = dotMatrix->operator()(j,kk);
225  x1 = output->operator[](kk)[0];
226  y1 = output->operator[](kk)[1];
227  z1 = output->operator[](kk)[2];
228  double value = x*x1 + y*y1 + z*z1;
229  if (value<0)
230  value = -value;
231  if (value > maxDots[j] )
232  maxDots[j] = value;
233  maxValue_j = maxDots[j];
234  }
235 
236  if (maxValue_j < minValue )
237  {
238  minValue = maxValue_j;
239  index = j;
240  }
241  }
242  }
243 
244  chosenIndex.push_back(index);
245  hasChosen[index] = 1;
246  indices->operator[](0).push_back(index);
247  }
248  // clock.Stop();
249  // std::cout << "after loop: " << clock.GetTotal() << std::endl;
250  }
251  else if (this->m_CriteriaType==Self::ELECTROSTATIC)
252  {
253  double x,y,z,x1,y1,z1;
254  bool orderEqual2 = std::abs(this->m_ElectrostaticOrder-2)<1e-8 ? true : false;
255  std::vector<double> energy(numberOfSamples,0.0);
256  // clock.Start();
257  for ( unsigned int i = (this->m_InitialOrientations?numberInitialSamples:1); i < totalNumbers; i += 1 )
258  {
259  int index=-1;
260  double minValue=std::numeric_limits<double>::max();
261  for ( unsigned int j = 0; j < numberOfSamples; j += 1 )
262  {
263  if (hasChosen[j]==0)
264  {
265  x = output->operator[](j)[0];
266  y = output->operator[](j)[1];
267  z = output->operator[](j)[2];
268  if (energy[j]<1e-10)
269  {
270  for ( unsigned int k = 0; k < chosenIndex.size(); k += 1 )
271  {
272  unsigned int kk = chosenIndex[k];
273  // double value = elecMatrix->operator()(j,kk);
274  x1 = output->operator[](kk)[0];
275  y1 = output->operator[](kk)[1];
276  z1 = output->operator[](kk)[2];
277  double dot = x*x1 + y*y1 + z*z1;
278  double result1 = 2+2*dot;
279  double result2 = 2-2*dot;
280  double value = (!orderEqual2) ? (1.0/std::pow(result1, this->m_ElectrostaticOrder*0.5) + 1.0/std::pow(result2, this->m_ElectrostaticOrder*0.5)) : (1.0/result1 + 1.0/result2);
281  energy[j] += value;
282  }
283  }
284  else
285  {
286  unsigned int kk = chosenIndex.back();
287  // double value = dotMatrix->operator()(j,kk);
288  x1 = output->operator[](kk)[0];
289  y1 = output->operator[](kk)[1];
290  z1 = output->operator[](kk)[2];
291  double dot = x*x1 + y*y1 + z*z1;
292  double result1 = 2+2*dot;
293  double result2 = 2-2*dot;
294  double value = (!orderEqual2) ? (1.0/std::pow(result1, this->m_ElectrostaticOrder*0.5) + 1.0/std::pow(result2, this->m_ElectrostaticOrder*0.5)) : (1.0/result1 + 1.0/result2);
295  energy[j] += value;
296  }
297  if (energy[j] < minValue )
298  {
299  minValue = energy[j];
300  index = j;
301  }
302  }
303  }
304  chosenIndex.push_back(index);
305  hasChosen[index] = 1;
306  indices->operator[](0).push_back(index);
307  }
308  // clock.Stop();
309  // std::cout << "after loop: " << clock.GetTotal() << std::endl;
310  }
311  else
312  utlGlobalException(true, "wrong m_CriteriaType");
313 
314  }
315  else
316  {
317 
318  if (this->m_CriteriaType==Self::DISTANCE)
319  {
320 
321  double x,y,z,x1,y1,z1;
322  std::vector< std::vector<double> > maxDotsVec;
323  for ( int i = 0; i < numberOfShells+1; i += 1 )
324  {
325  std::vector<double> dotTemp(numberOfSamples,-1);
326  maxDotsVec.push_back(dotTemp);
327  }
328 
329  /**************************************************************************************/
330  // first several points for the shell without initialization.
331  unsigned startI = (this->m_InitialOrientations?numberInitialSamples:1);
332  for ( ; startI < numberOfShells; startI += 1 )
333  {
334  int index=-1;
335  double minValue=std::numeric_limits<double>::max();
336  for ( unsigned int j = 0; j < numberOfSamples; j += 1 )
337  {
338  if (hasChosen[j]==0)
339  {
340  double maxValue_j = -std::numeric_limits<double>::max();
341  x = output->operator[](j)[0];
342  y = output->operator[](j)[1];
343  z = output->operator[](j)[2];
344  if (maxDotsVec[numberOfShells][j]< 0 )
345  {
346  for ( unsigned int k = 0; k < chosenIndex.size(); k += 1 )
347  {
348  unsigned int kk = chosenIndex[k];
349  // double value = dotMatrix->operator()(j,kk);
350  x1 = output->operator[](kk)[0];
351  y1 = output->operator[](kk)[1];
352  z1 = output->operator[](kk)[2];
353  double value = x*x1 + y*y1 + z*z1;
354  if (value<0)
355  value = -value;
356  if (value>maxValue_j)
357  maxValue_j = value;
358  }
359  maxDotsVec[numberOfShells][j] = maxValue_j;
360  }
361  else
362  {
363  unsigned int kk = chosenIndex.back();
364  // double value = dotMatrix->operator()(j,kk);
365  x1 = output->operator[](kk)[0];
366  y1 = output->operator[](kk)[1];
367  z1 = output->operator[](kk)[2];
368  double value = x*x1 + y*y1 + z*z1;
369  if (value<0)
370  value = -value;
371  if (value > maxDotsVec[numberOfShells][j] )
372  maxDotsVec[numberOfShells][j] = value;
373  maxValue_j = maxDotsVec[numberOfShells][j];
374  }
375 
376  if (maxValue_j < minValue )
377  {
378  minValue = maxValue_j;
379  index = j;
380  }
381  }
382  }
383  chosenIndex.push_back(index);
384  hasChosen[index] = 1;
385  IndexVectorType tmp;
386  tmp.push_back(index);
387  indices->push_back(tmp);
388  }
389 
390  utlException(indices->size()!=numberOfShells, "Logic ERROR!");
391  // utlPrintVar1(true, startI);
392  // utl::Save2DVector(*indices, std::cout<<"indices = ");
393 
394  /**************************************************************************************/
395  // other points
396  int chosedShellIndex_final = -1;
397  int chosedShellIndex_final_last = -1;
398  for ( unsigned int i = startI; i < totalNumbers; i += 1 )
399  {
400  int index=-1;
401  double maxDist = -std::numeric_limits<double>::max();
402  int chosedShellIndex = -1;
403 
404  for ( unsigned int j = 0; j < numberOfSamples; j += 1 )
405  {
406  if (hasChosen[j]==0)
407  {
408 
409  // maximal dots in single shells
410  std::vector<double> maxValue_j(numberOfShells+1, -std::numeric_limits<double>::max());
411  x = output->operator[](j)[0];
412  y = output->operator[](j)[1];
413  z = output->operator[](j)[2];
414 
415  // utlPrintVar1(true, chosedShellIndex_final);
416  if (chosedShellIndex_final_last>=0)
417  {
418  // update stored maxDotsVec based on last added point
419  unsigned int kk = chosenIndex.back();
420  // double value = dotMatrix->operator()(j,kk);
421  x1 = output->operator[](kk)[0];
422  y1 = output->operator[](kk)[1];
423  z1 = output->operator[](kk)[2];
424  double value = x*x1 + y*y1 + z*z1;
425  if (value<0)
426  value = -value;
427  if (value > maxDotsVec[chosedShellIndex_final_last][j] )
428  maxDotsVec[chosedShellIndex_final_last][j] = value;
429  if (value > maxDotsVec[numberOfShells][j] )
430  maxDotsVec[numberOfShells][j] = value;
431  }
432 
433  for ( unsigned int s = 0; s < numberOfShells; s += 1 )
434  {
435  int numbers = (*indices)[s].size();
436  if (numbers < this->m_NumbersInShell[s])
437  {
438  // only add point if there is no enough number
439  if (maxDotsVec[s][j]<0)
440  {
441  // only calculate all dots in the first begining
442  for ( unsigned int k = 0; k < numbers; k += 1 )
443  {
444  unsigned int kk = (*indices)[s][k];
445  // double value = dotMatrix->operator()(j,kk);
446  x1 = output->operator[](kk)[0];
447  y1 = output->operator[](kk)[1];
448  z1 = output->operator[](kk)[2];
449  double value = x*x1 + y*y1 + z*z1;
450  if (value<0)
451  value = -value;
452  if (value>maxValue_j[s])
453  maxValue_j[s] = value;
454  }
455  maxDotsVec[s][j] = maxValue_j[s];
456  if (maxValue_j[s]>maxDotsVec[numberOfShells][j])
457  maxDotsVec[numberOfShells][j] = maxValue_j[s];
458  }
459  else
460  {
461  maxValue_j[s] = maxDotsVec[s][j];
462  }
463  }
464  }
465 
466  // maximal dot in whole shell, which is the maxinmal value amoung maximal dots in each shell
467  // for ( unsigned int s = 0; s < numberOfShells; s += 1 )
468  // {
469  // if (maxDotsVec[s][j] > maxValue_j[numberOfShells])
470  // maxValue_j[numberOfShells] = maxDotsVec[s][j];
471  // }
472 
473  // maximal dot in whole shell
474  maxValue_j[numberOfShells] = maxDotsVec[numberOfShells][j];
475 
476 
477  // choose the shell with minimal maximal dot
478  chosedShellIndex = -1;
479  double chosedShellDot = std::numeric_limits<double>::max();
480  for ( unsigned int s = 0; s < numberOfShells; s += 1 )
481  {
482  if ( (*indices)[s].size()==this->m_NumbersInShell[s])
483  continue;
484  if (maxValue_j[s] < chosedShellDot)
485  {
486  chosedShellIndex = s;
487  chosedShellDot = maxValue_j[s];
488  }
489  }
490  // utlPrintVar2(true, j, chosedShellIndex);
491  // utl::PrintVector(maxValue_j, "maxValue_j");
492 
493  // distance in cost function which combines whole shell and the chosed shell
494  double costDist = (1-this->m_WeightForSingleShell)*std::acos(maxValue_j[numberOfShells])
495  + this->m_WeightForSingleShell*std::acos(maxValue_j[chosedShellIndex]);
496  if (costDist > maxDist )
497  {
498  maxDist = costDist;
499  index = j;
500  chosedShellIndex_final = chosedShellIndex;
501  }
502 
503  }
504  }
505 
506  chosenIndex.push_back(index);
507  hasChosen[index] = 1;
508  indices->operator[](chosedShellIndex_final).push_back(index);
509  chosedShellIndex_final_last = chosedShellIndex_final;
510  }
511 
512  }
513  else if (this->m_CriteriaType==Self::ELECTROSTATIC)
514  {
515 
516  utlGlobalException(true, "TODO");
517  }
518  else
519  utlGlobalException(true, "wrong m_CriteriaType");
520 
521  }
522 }
523 
524 }
525 
526 
527 #endif
528 
529 
530 
utl.h
helper functions specifically used in dmritool
utl::ConnectUtlMatrix
NDArray< T, 2 > ConnectUtlMatrix(const NDArray< T, 2 > &m1, const NDArray< T, 2 > &m2, const bool isConnectRow)
Definition: utlNDArrayFunctions.h:1033
utlException
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
utl::max
const T & max(const T &a, const T &b)
Return the maximum between a and b.
Definition: utlCore.h:263
spams::abs
T abs(const T x)
template version of the fabs function
utlGlobalException
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
CARTESIAN_TO_CARTESIAN
Definition: utlDMRI.h:21
DIRECTION_NODUPLICATE
Definition: utlDMRI.h:35
Generated by   doxygen 1.8.11