DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkDWIReader.hxx
Go to the documentation of this file.
1 
19 #ifndef __itkDWIReader_hxx
20 #define __itkDWIReader_hxx
21 
22 #include "itkDWIReader.h"
23 #include "itkImageFileReader.h"
24 #include "itkImageRegionIteratorWithIndex.h"
25 #include "utl.h"
26 
27 
28 namespace itk
29 {
30 
31 template <class TPixelType, unsigned int VImageDimension>
32 DWIReader<TPixelType, VImageDimension>
33 ::DWIReader()
34 {
35  m_NormalizeDWI = true;
36  // m_DWIWithB0 = false;
37  m_IsInput4DImage = true;
38  m_CorrectDWIValues = true;
39  m_ShowWarnings = true;
40 
41  m_MaskImage = MaskImageType::New();
42  m_B0Image = B0ImageType::New();
43 
44  m_SamplingSchemeQSpace = SamplingSchemeQSpaceType::New();
45 }
46 
47 template <class TPixelType, unsigned int VImageDimension>
48 bool
49 DWIReader<TPixelType, VImageDimension>
50 ::DetermineIsInput4DImage(const std::string& dataStr)
51 {
52  typedef Image<TPixelType,4> DWI4DImageType;
53  typename DWI4DImageType::Pointer dwi4DTempImage;
54  itk::ReadImageInformation<DWI4DImageType>(dataStr, dwi4DTempImage);
55  typename DWI4DImageType::SizeType size = dwi4DTempImage->GetLargestPossibleRegion().GetSize();
56  int numberOfDWI = size[3];
57  if (numberOfDWI>1)
58  return true; // 4D Image
59  else
60  {
61  typename DWIImageType::Pointer dwiTempImage;
62  itk::ReadImageInformation<DWIImageType>(dataStr, dwiTempImage);
63  numberOfDWI = dwiTempImage->GetNumberOfComponentsPerPixel();
64  if (numberOfDWI>1)
65  return false; // VectorImage
66  else
67  return true; // 4D Image, one component
68  }
69 }
70 
71 template <class TPixelType, unsigned int VImageDimension>
72 void
73 DWIReader<TPixelType, VImageDimension>
74 ::ReadFromConfigurationFile(const std::string& file)
75 {
76  utlShowPosition(this->GetDebug());
77 
78  // read file
79  utlGlobalException(m_ConfigurationFile=="", "need to set the configuration file!");
80  std::vector<std::vector<std::string> > stringMatrix;
81  utl::ReadLinesFirstlineCheck(m_ConfigurationFile, stringMatrix);
82 
83  // check configuration file and dwi size
84  std::string dwiPath, dwiFile;
85  utl::GetPath(file, dwiPath, dwiFile);
86  double bvalue=-1;
87  std::string dataStr, gradStr, indexStr, bStr;
88  std::vector<std::vector<std::string> > gradStrTempMatrix;
89  std::vector<int> numDWINob0, additionalB0;
90  std::vector<std::vector<double> > bVec;
91  std::vector<std::vector<std::vector<std::string> > > gradStrMatrix;
92  std::vector<std::vector<int> > b0IndexVec;
93  std::vector<std::vector<int> > selectIndexVec;
94 
95  typedef ImageFileReader<DWIImageType> DWIReaderType;
96  typename DWIImageType::Pointer dwiTempImage;
97 
98  typedef Image<TPixelType, 4> DWI4DImageType;
99  typedef ImageFileReader<DWI4DImageType> DWI4DReaderType;
100  typename DWI4DImageType::Pointer dwi4DTempImage;
101 
102  bool isB0Set = !IsImageEmpty(m_B0Image);
103  double bThresholdSingleShell = m_SamplingSchemeQSpace->GetBThresholdSingleShell();
104 
105  // Read Input File Info
106  int numberOfDWIsWithoutB0=0;
107  int numberOfB0=0;
108  std::vector<std::string> ossStr;
109  std::string dataStr0;
110  for ( int i = 0; i < stringMatrix.size(); i += 1 )
111  {
112  // utlPrintVar1(true, i);
113  utlGlobalException(stringMatrix[i].size()!=3 && stringMatrix[i].size()!=4, "wrong size of configuration file");
114 
115  dataStr = utl::IsFileExist(stringMatrix[i][2]) ? stringMatrix[i][2] : dwiPath+stringMatrix[i][2] ;
116  gradStr = utl::IsFileExist(stringMatrix[i][1]) ? stringMatrix[i][1] : dwiPath+stringMatrix[i][1] ;
117 
118  if (i==0)
119  {
120  dataStr0 = utl::IsFileExist(stringMatrix[i][2]) ? stringMatrix[i][2] : dwiPath+stringMatrix[i][2];
121  utlGlobalException( !IsImageEmpty(m_MaskImage) && !itk::VerifyImageSize<MaskImageType>(m_MaskImage, dataStr0, true), "wrong size of m_MaskImage and DWI file " << stringMatrix[i][2] << ". Inconsistent information!");
122  utlGlobalException( !IsImageEmpty(m_B0Image) && !itk::VerifyImageSize<B0ImageType>(m_B0Image, dataStr0, true), "wrong size of b0Image and DWI file " << stringMatrix[i][2] << ". Inconsistent information! m_B0Image="<<m_B0Image);
123  m_IsInput4DImage = DetermineIsInput4DImage(dataStr0);
124  }
125  else
126  {
127  utlGlobalException( !itk::VerifyImageInformation(dataStr0, dataStr, true), "DWI files "<< dataStr << " and " << dataStr0 << " have inconsistent information!");
128  }
129 
130  if (stringMatrix[i].size()==4)
131  {
132  indexStr = utl::IsFileExist(stringMatrix[i][3]) ? stringMatrix[i][3] : dwiPath+stringMatrix[i][3];
133  std::vector<int> tmpVec;
134  utl::ReadVector(indexStr,tmpVec);
135  selectIndexVec.push_back(tmpVec);
136  }
137 
138  utl::ReadLinesFirstlineCheck(gradStr, gradStrTempMatrix);
139  std::vector<double> bVecTemp;
140  if (utl::IsNumber(stringMatrix[i][0]))
141  {
142  std::istringstream ( stringMatrix[i][0] ) >> bvalue;
143  bVecTemp.clear();
144  bVecTemp.assign(gradStrTempMatrix.size(), bvalue);
145  }
146  else
147  {
148  bStr = utl::IsFileExist(stringMatrix[i][0]) ? stringMatrix[i][0] : dwiPath+stringMatrix[i][0];
149  utl::ReadVector(bStr,bVecTemp);
150  utlGlobalException(gradStrTempMatrix.size()!=gradStrTempMatrix.size(), "wrong size of b values in " << bStr << " and gradients in " << gradStr);
151  }
152  bVec.push_back(bVecTemp);
153 
154  std::vector<int> b0IndexTemp;
155  for ( int j = 0; j < gradStrTempMatrix.size(); j += 1 )
156  {
157  utlGlobalException(gradStrTempMatrix[j].size()!=3, "wrong size of gradients " << gradStr);
158  if (stringMatrix[i].size()==4 && !utl::IsInVector(selectIndexVec[i],j))
159  continue;
160  double gradx = utl::ConvertStringToNumber<double>(gradStrTempMatrix[j][0]);
161  double grady = utl::ConvertStringToNumber<double>(gradStrTempMatrix[j][1]);
162  double gradz = utl::ConvertStringToNumber<double>(gradStrTempMatrix[j][2]);
163  bool gradIsZero = std::abs(gradx)<1e-4 && std::abs(grady)<1e-4 && std::abs(gradz)<1e-4;
164  bool bIsZero = std::abs(bVecTemp[j])<(bThresholdSingleShell>0?bThresholdSingleShell:1e-4);
165  // std::cout << "bThresholdSingleShell=" << bThresholdSingleShell << ", bIsZero="<<bIsZero << std::endl << std::flush;
166  // utlGlobalException(gradIsZero && !bIsZero || !gradIsZero && bIsZero, "Wrong logic in gradients and b values! bValue = "<< bVecTemp[j] << ", grad=" <<"("<<gradx<<","<<grady<<","<<gradz<<").");
167  if (gradIsZero || bIsZero)
168  b0IndexTemp.push_back(j);
169  }
170  utlGlobalException(isB0Set && b0IndexTemp.size()>0, "Since b0 is mannually set already, the configuration file should not have b0 image");
171  gradStrMatrix.push_back(gradStrTempMatrix);
172  b0IndexVec.push_back(b0IndexTemp);
173 
174  std::ostringstream oss;
175  if (utl::IsNumber(stringMatrix[i][0]))
176  oss << "> b="<< bvalue;
177  else
178  oss << "> bStr="<< bStr;
179  oss << ", gradStrTempMatrix.size()=" << gradStrTempMatrix.size() << ", 4D image = " << dataStr;
180  if (stringMatrix[i].size()==4)
181  oss << ", index=" << indexStr;
182  // oss << std::endl << std::flush;
183  ossStr.push_back(oss.str());
184 
185 
186  int numberOfDWI;
187  if (m_IsInput4DImage)
188  {
189  itk::ReadImageInformation<DWI4DImageType>(dataStr, dwi4DTempImage);
190  typename DWI4DImageType::SizeType size = dwi4DTempImage->GetLargestPossibleRegion().GetSize();
191  numberOfDWI = size[3];
192  }
193  else
194  {
195  itk::ReadImageInformation<DWIImageType>(dataStr, dwiTempImage);
196  numberOfDWI = dwiTempImage->GetNumberOfComponentsPerPixel();
197  }
198 
199  utlGlobalException(numberOfDWI<gradStrTempMatrix.size(), "The number of gradients in " << gradStr << " is "<< gradStrTempMatrix.size() << " more than number of DWIs in " << dataStr << "="<<numberOfDWI);
200  additionalB0.push_back(numberOfDWI - gradStrTempMatrix.size());
201  if (stringMatrix[i].size()==3)
202  {
203  numDWINob0.push_back(gradStrTempMatrix.size() - b0IndexTemp.size());
204  numberOfDWIsWithoutB0 += gradStrTempMatrix.size() - b0IndexTemp.size();
205  numberOfB0 += b0IndexTemp.size() + numberOfDWI - gradStrTempMatrix.size();
206  }
207  else
208  {
209  numDWINob0.push_back(selectIndexVec[i].size() - b0IndexTemp.size());
210  numberOfDWIsWithoutB0 += selectIndexVec[i].size() - b0IndexTemp.size();
211  numberOfB0 += b0IndexTemp.size() + numberOfDWI - gradStrTempMatrix.size();
212  }
213 
214  // numDWInob0 records the start index for i-th dwi file
215  if (i>0)
216  numDWINob0[i] += numDWINob0[i-1];
217  }
218 
219  // utlGlobalException( !this->GetB0Image() && numberOfB0==0, "no b0 images in " << file);
220  std::cout << numberOfB0 << " b0 images, " << numberOfDWIsWithoutB0 << " DWIs without b0. "<< std::endl << std::flush;
221  if (this->GetDebug())
222  {
223  utl::PrintVector(numDWINob0, "numDWINob0");
224  utl::PrintVector(additionalB0, "additionalB0");
225  for ( int jj = 0; jj < b0IndexVec.size(); jj += 1 )
226  utl::PrintVector(b0IndexVec[jj], "b0IndexVec[jj]");
227  }
228 
229  // Read data in input file
230  typename DWIImageType::Pointer dwiImage = this->GetOutput();
231  STDVectorPointer bVector = m_SamplingSchemeQSpace->GetBVector();
232  MatrixPointer orientationsCartesian = m_SamplingSchemeQSpace->GetOrientationsCartesian();
233  bVector->clear();
234  orientationsCartesian->ReSize(numberOfDWIsWithoutB0, 3);
235  typename B0ImageType::IndexType b0Index;
236  typename B0ImageType::PixelType b0Pixel;
237  PixelType dwiPixel, dwiZeroPixel;
238  dwiPixel.SetSize(numberOfDWIsWithoutB0);
239  dwiZeroPixel.SetSize(numberOfDWIsWithoutB0);
240  dwiZeroPixel.Fill(0.0);
241 
242  for ( int i = 0; i < stringMatrix.size(); i += 1 )
243  {
244  if (stringMatrix.size()==1)
245  std::cout << "loading the " << i+1 << "th DWI data from total " << stringMatrix.size() << " DWI data" << std::endl;
246  std::cout << ossStr[i] << std::endl << std::flush;
247  dataStr = utl::IsFileExist(stringMatrix[i][2]) ? stringMatrix[i][2] : dwiPath+stringMatrix[i][2];
248 
249  // b vector, gradients
250  // utl::PrintVector(bVec[i], "bVec[i]");
251  for ( int j = 0; j < bVec[i].size(); j += 1 )
252  {
253  if (!utl::IsInVector(b0IndexVec[i], j) && ( stringMatrix[i].size()==3 || (stringMatrix[i].size()==4 && utl::IsInVector(selectIndexVec[i],j))) )
254  {
255  // utlPrintVar1(true, j);
256  bVector->push_back(bVec[i][j]);
257  double gx = utl::ConvertStringToNumber<double>(gradStrMatrix[i][j][0]);
258  double gy = utl::ConvertStringToNumber<double>(gradStrMatrix[i][j][1]);
259  double gz = utl::ConvertStringToNumber<double>(gradStrMatrix[i][j][2]);
260  double gnorm = std::sqrt(gx*gx+gy*gy+gz*gz);
261  (*orientationsCartesian)(bVector->size()-1,0) = gx/gnorm;
262  (*orientationsCartesian)(bVector->size()-1,1) = gy/gnorm;
263  (*orientationsCartesian)(bVector->size()-1,2) = gz/gnorm;
264  // utlPrintVar3(true, (*orientationsCartesian)(bVector->size()-1,0), (*orientationsCartesian)(bVector->size()-1,1), (*orientationsCartesian)(bVector->size()-1,2));
265  }
266  }
267 
268  // read individual dwi
269  int numberOfDWI;
270  if (m_IsInput4DImage)
271  {
272  itk::ReadImage<DWI4DImageType>(dataStr, dwi4DTempImage);
273  typename DWI4DImageType::SizeType size = dwi4DTempImage->GetLargestPossibleRegion().GetSize();
274  numberOfDWI = size[3];
275  }
276  else
277  {
278  itk::ReadImage<DWIImageType>(dataStr, dwiTempImage);
279  numberOfDWI = dwiTempImage->GetNumberOfComponentsPerPixel();
280  }
281 
282  // allocate b0, dwi
283  if (i==0)
284  {
285  if (IsImageEmpty(m_B0Image))
286  {
287  m_B0Image = B0ImageType::New();
288  if (m_IsInput4DImage)
289  itk::CopyImageInformation<DWI4DImageType, B0ImageType>( dwi4DTempImage, m_B0Image);
290  else
291  itk::CopyImageInformation<DWIImageType, B0ImageType>( dwiTempImage, m_B0Image);
292  m_B0Image->Allocate();
293  if (numberOfB0>0)
294  m_B0Image->FillBuffer(0);
295  else
296  {
297  std::cout << "b0=1 is used" << std::endl << std::flush;
298  m_B0Image->FillBuffer(1.0); // assume b0=1.0, if there is no b0 images.
299  }
300  }
301  if (m_IsInput4DImage)
302  itk::CopyImageInformation<DWI4DImageType, DWIImageType>( dwi4DTempImage, dwiImage);
303  else
304  itk::CopyImageInformation<DWIImageType, DWIImageType>( dwiTempImage, dwiImage);
305  dwiImage->SetNumberOfComponentsPerPixel( numberOfDWIsWithoutB0 );
306  dwiImage->Allocate();
307  dwiImage->FillBuffer(dwiZeroPixel);
308  }
309 
310  // concatenate dwi
311  ImageRegionIteratorWithIndex<B0ImageType> b0It(m_B0Image, m_B0Image->GetLargestPossibleRegion());
312  ImageRegionIteratorWithIndex<DWIImageType> dwiIt(dwiImage, dwiImage->GetLargestPossibleRegion());
313  ImageRegionIteratorWithIndex<MaskImageType> maskIt;
314  if (!IsImageEmpty(m_MaskImage))
315  maskIt = ImageRegionIteratorWithIndex<MaskImageType>(m_MaskImage, m_MaskImage->GetLargestPossibleRegion());
316  IndexType dwiIndex;
317 
318  if (m_IsInput4DImage)
319  {
320  typename DWI4DImageType::IndexType dwi4DTempIndex;
321  for (b0It.GoToBegin(), dwiIt.GoToBegin(), maskIt.GoToBegin();
322  !b0It.IsAtEnd();
323  ++b0It, ++dwiIt, ++maskIt)
324  {
325  if (!IsImageEmpty(m_MaskImage) && maskIt.Get()<=1e-10)
326  {
327  dwiIt.Set(dwiZeroPixel);
328  b0It.Set(0.0);
329  continue;
330  }
331 
332  b0Index = b0It.GetIndex();
333  if (this->GetDebug())
334  std::cout << "index = " << b0Index << std::endl << std::flush;
335  for ( int kk = 0; kk < VImageDimension; kk += 1 )
336  dwi4DTempIndex[kk] = b0Index[kk];
337 
338  std::ostringstream ossWarn, ossWarnTmp;
339  // additional b0
340  double sumPixelb=0;
341  for ( int kk = 0; kk < additionalB0[i]; kk += 1 )
342  {
343  dwi4DTempIndex[VImageDimension] = kk;
344  double dwi4DTempPixel = dwi4DTempImage->GetPixel(dwi4DTempIndex);
345  sumPixelb += dwi4DTempPixel*dwi4DTempPixel;
346  if (dwi4DTempPixel<1e-10)
347  ossWarnTmp << "zero b0 at "<< dwi4DTempIndex <<", " << "b0=" << dwi4DTempPixel << ", ";
348  if (numberOfB0>0)
349  {
350  b0Pixel = b0It.Get() + dwi4DTempPixel/((double)numberOfB0);
351  b0It.Set(b0Pixel);
352  }
353  }
354 
355  // b0 in gradients
356  for ( int kk = 0; kk < b0IndexVec[i].size(); kk += 1 )
357  {
358  dwi4DTempIndex[VImageDimension] = b0IndexVec[i][kk]+additionalB0[i];
359  double dwi4DTempPixel = dwi4DTempImage->GetPixel(dwi4DTempIndex);
360  sumPixelb += dwi4DTempPixel*dwi4DTempPixel;
361  if (dwi4DTempPixel<1e-10)
362  ossWarnTmp << "zero b0 at "<< dwi4DTempIndex <<", " << "b0=" << dwi4DTempPixel << ", ";
363  if (numberOfB0>0)
364  {
365  b0Pixel = b0It.Get() + dwi4DTempPixel/((double)numberOfB0);
366  b0It.Set(b0Pixel);
367  }
368  }
369 
370  if (b0It.Get()<1e-10)
371  {
372  dwiIt.Set(dwiZeroPixel);
373  continue;
374  }
375 
376  // dwi without b0
377  double sumPixeldwi=0;
378  dwiPixel = dwiIt.Get();
379  int jj=0;
380  for ( int kk = additionalB0[i]; kk < numberOfDWI; kk += 1 )
381  {
382  if (!utl::IsInVector(b0IndexVec[i], kk-additionalB0[i]) && ( stringMatrix[i].size()==3 || (stringMatrix[i].size()==4 && utl::IsInVector(selectIndexVec[i],kk-additionalB0[i]))) )
383  {
384  dwi4DTempIndex[VImageDimension] = kk;
385  double dwi4DTempPixel = dwi4DTempImage->GetPixel(dwi4DTempIndex);
386  sumPixeldwi += dwi4DTempPixel*dwi4DTempPixel;
387  if (dwi4DTempPixel<1e-10)
388  ossWarnTmp << "zero dwi at "<< dwi4DTempIndex <<", ";
389  dwiPixel[jj + (i==0?0:numDWINob0[i-1])] = dwi4DTempPixel;
390  jj++;
391  }
392  }
393  if (ossWarnTmp.str()!="" && sumPixelb>0 && sumPixeldwi>0)
394  ossWarn << ossWarnTmp.str();
395 
396  if (m_ShowWarnings && (b0It.Get()>1e-10 || dwiPixel.GetSquaredNorm()>1e-10) && ossWarn.str()!="")
397  {
398  // utlWarning(true, "Warning: " << ossWarn.str() << "dwiPixel=" << dwiPixel<<";");
399  std::cout << "Warning: " << ossWarn.str() << std::endl << std::flush;
400  }
401 
402  if (this->GetDebug())
403  {
404  utlPrintVar(true, b0It.Get());
405  itk::PrintVariableLengthVector(dwiPixel, "dwiPixel");
406  }
407 
408  dwiIt.Set(dwiPixel);
409 
410  }
411 
412  }
413  else
414  {
415  ImageRegionIteratorWithIndex<DWIImageType> dwiTempIt(dwiTempImage, dwiTempImage->GetLargestPossibleRegion());
416  PixelType dwiTempPixel;
417  for (b0It.GoToBegin(), dwiIt.GoToBegin(), dwiTempIt.GoToBegin(), maskIt.GoToBegin();
418  !b0It.IsAtEnd();
419  ++b0It, ++dwiIt, ++dwiTempIt, ++maskIt)
420  {
421  if (!IsImageEmpty(m_MaskImage) && maskIt.Get()<=1e-10)
422  {
423  dwiIt.Set(dwiZeroPixel);
424  b0It.Set(0.0);
425  continue;
426  }
427 
428  b0Index = b0It.GetIndex();
429  if (this->GetDebug())
430  std::cout << "index = " << b0Index << std::endl << std::flush;
431 
432  dwiTempPixel = dwiTempIt.Get();
433  std::ostringstream ossWarn;
434  // additional b0
435  for ( int kk = 0; kk < additionalB0[i]; kk += 1 )
436  {
437  if (dwiTempPixel[kk]<1e-10)
438  ossWarn << "zero b0 at ["<< b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<","<<kk <<"], b0=" << dwiTempPixel[kk] << "\n";
439  if (numberOfB0>0)
440  {
441  b0Pixel = b0It.Get() + dwiTempPixel[kk]/((double)numberOfB0);
442  b0It.Set(b0Pixel);
443  }
444  }
445  // b0 in gradients
446  for ( int kk = 0; kk < b0IndexVec[i].size(); kk += 1 )
447  {
448  if (dwiTempPixel[kk]<1e-10)
449  ossWarn << "zero b0 at ["<< b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<","<<kk <<"], b0=" << dwiTempPixel[kk] << "\n";
450  if (numberOfB0>0)
451  {
452  b0Pixel = b0It.Get() + dwiTempPixel[ b0IndexVec[i][kk]+additionalB0[i] ]/((double)numberOfB0);
453  b0It.Set(b0Pixel);
454  }
455  }
456 
457  if (b0It.Get()<1e-10)
458  {
459  dwiIt.Set(dwiZeroPixel);
460  continue;
461  }
462 
463  // dwi without b0
464  dwiPixel = dwiIt.Get();
465  int jj=0;
466  for ( int kk = additionalB0[i]; kk < numberOfDWI; kk += 1 )
467  {
468  if (!utl::IsInVector(b0IndexVec[i], kk-additionalB0[i]) && ( stringMatrix[i].size()==3 || (stringMatrix[i].size()==4 && utl::IsInVector(selectIndexVec[i],kk-additionalB0[i]))) )
469  {
470  if (dwiTempPixel[kk]<1e-10)
471  ossWarn << "zero dwi at ["<< b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<","<<kk <<"], dwi value=" << dwiTempPixel[kk] << "\n";
472  dwiPixel[jj + (i==0?0:numDWINob0[i-1])] = dwiTempPixel[kk];
473  jj++;
474  }
475  }
476 
477  if (m_ShowWarnings && (b0It.Get()>1e-10 && dwiPixel.GetSquaredNorm()>1e-10) && ossWarn.str()!="")
478  {
479  // utlWarning(true, "Warning: " << ossWarn.str() << "dwiPixel=" << dwiPixel<<";");
480  std::cout << "Warning: " << ossWarn.str() << std::endl << std::flush;
481  }
482 
483  if (this->GetDebug())
484  {
485  utlPrintVar(true, b0It.Get());
486  itk::PrintVariableLengthVector(dwiPixel, "dwiPixel");
487  }
488  dwiIt.Set(dwiPixel);
489  }
490 
491  }
492 
493  }
494 
495  m_SamplingSchemeQSpace->SetOrientationsCartesian(orientationsCartesian);
496  m_SamplingSchemeQSpace->SetBVector(bVector);
497 }
498 
499 template <class TPixelType, unsigned int VImageDimension>
500 void
501 DWIReader<TPixelType, VImageDimension>
502 ::CorrectDWI()
503 {
504  utlShowPosition(this->GetDebug());
505  utlGlobalException(!m_B0Image, "no b0 images");
506  typename DWIImageType::Pointer dwiImage = this->GetOutput();
507 
508  utlGlobalException(!dwiImage, "no dwi images");
509  utlGlobalException(!(itk::VerifyImageSize<DWIImageType, B0ImageType>(dwiImage, m_B0Image, true)), "dwi image and b0 image have different information");
510 
511  ImageRegionIteratorWithIndex<B0ImageType> b0It(m_B0Image, m_B0Image->GetLargestPossibleRegion());
512  ImageRegionIteratorWithIndex<DWIImageType> dwiIt(dwiImage, dwiImage->GetLargestPossibleRegion());
513  ImageRegionIteratorWithIndex<MaskImageType> maskIt;
514  if (!IsImageEmpty(m_MaskImage))
515  maskIt = ImageRegionIteratorWithIndex<MaskImageType>(this->m_MaskImage, m_MaskImage->GetLargestPossibleRegion());
516  typename B0ImageType::IndexType b0Index;
517  typename B0ImageType::PixelType b0Pixel;
518  PixelType dwiPixel;
519  int numberOfDWIsWithoutB0 = dwiImage->GetNumberOfComponentsPerPixel();
520  dwiPixel.SetSize(numberOfDWIsWithoutB0);
521  for (b0It.GoToBegin(), dwiIt.GoToBegin(), maskIt.GoToBegin();
522  !b0It.IsAtEnd();
523  ++b0It, ++dwiIt, ++maskIt)
524  {
525  if (!IsImageEmpty(m_MaskImage) && maskIt.Get()<=1e-10)
526  continue;
527  b0Pixel = b0It.Get();
528  if (b0Pixel<1e-10)
529  {
530  dwiPixel = dwiIt.Get();
531  dwiPixel.Fill(0.0);
532  dwiIt.Set(dwiPixel);
533  continue;
534  }
535 
536  dwiPixel = dwiIt.Get();
537  if (dwiPixel.GetSquaredNorm()<1e-10)
538  {
539  b0It.Set(0.0);
540  continue;
541  }
542 
543  b0Index = b0It.GetIndex();
544  double maxData=0.0, minData=std::numeric_limits<double>::max();
545  for ( int j = 0; j < numberOfDWIsWithoutB0; j += 1 )
546  {
547  if (dwiPixel[j]<=0.99*b0Pixel && dwiPixel[j]>=1e-8*b0Pixel && dwiPixel[j]>maxData)
548  maxData = dwiPixel[j];
549  if (dwiPixel[j]<=0.99*b0Pixel && dwiPixel[j]>=1e-8*b0Pixel && dwiPixel[j]<minData)
550  minData = dwiPixel[j];
551  }
552 
553  if (maxData < minData)
554  {
555  if (this->GetDebug())
556  {
557  std::cout << "Logical ERROR: ("<<b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<"), all dwi values are not in [0, 0.99]*b0. " << std::endl << std::flush;
558  b0It.Set(0.0);
559  dwiPixel = dwiIt.Get();
560  dwiPixel.Fill(0.0);
561  dwiIt.Set(dwiPixel);
562  continue;
563  }
564  }
565 
566  for ( int j = 0; j < numberOfDWIsWithoutB0; j += 1 )
567  {
568  if (dwiPixel[j]>=(1+1e-8)*b0Pixel)
569  {
570  if (this->GetDebug())
571  {
572  std::cout << "Logical ERROR in itk::DWIReader, dwi("<<b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<","<<j<<")=" << dwiPixel[j] << " > " << "b0("<<b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<")=" << b0Pixel
573  << " Thus we force data["<<j<<"]="<< maxData <<", which is the maximal value of the other dwi values."<< std::endl;
574  }
575  dwiPixel[j] = maxData;
576  }
577  if (dwiPixel[j]<=1e-8*b0Pixel)
578  {
579  if (this->GetDebug())
580  {
581  std::cout << "Logical ERROR in itk::DWIReader, dwi("<<b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<","<<j<<")=" << dwiPixel[j] << " < " << "b0("<<b0Index[0]<<","<<b0Index[1]<<","<<b0Index[2]<<")=" << b0Pixel
582  << " Thus we force data["<<j<<"]="<< minData <<", which is the minimal value of the other dwi values."<< std::endl;
583  }
584  dwiPixel[j] = minData;
585  }
586  }
587  dwiIt.Set(dwiPixel);
588  }
589 
590  this->Modified();
591 }
592 
593 template <class TPixelType, unsigned int VImageDimension>
594 void
595 DWIReader<TPixelType, VImageDimension>
596 ::NormalizeDWI()
597 {
598  utlShowPosition(this->GetDebug());
599  utlGlobalException(!m_B0Image, "no b0 images");
600  typename DWIImageType::Pointer dwiImage = this->GetOutput();
601 
602  utlGlobalException(!dwiImage, "no dwi images");
603  utlGlobalException(!(itk::VerifyImageSize<DWIImageType, B0ImageType>(dwiImage, m_B0Image, true)), "dwi image and b0 image have different information");
604 
605  ImageRegionIterator<B0ImageType> b0It(m_B0Image, m_B0Image->GetLargestPossibleRegion());
606  ImageRegionIterator<DWIImageType> dwiIt(dwiImage, dwiImage->GetLargestPossibleRegion());
607  typename B0ImageType::PixelType b0Pixel;
608  PixelType dwiPixel;
609  dwiPixel.SetSize(dwiImage->GetNumberOfComponentsPerPixel());
610 
611  for (b0It.GoToBegin(), dwiIt.GoToBegin();
612  !b0It.IsAtEnd();
613  ++b0It, ++dwiIt)
614  {
615  b0Pixel = b0It.Get();
616  if (b0Pixel<1e-10)
617  continue;
618  dwiPixel = dwiIt.Get();
619  dwiPixel /= b0Pixel;
620  dwiIt.Set(dwiPixel);
621  }
622  this->Modified();
623 }
624 
625 template <class TPixelType, unsigned int VImageDimension>
626 void
627 DWIReader<TPixelType, VImageDimension>
628 ::GenerateData()
629 {
630  if (m_ConfigurationFile!="")
631  ReadFromConfigurationFile(m_ConfigurationFile);
632  else
633  utlGlobalException(true, "unknown format! ");
634 
635  if (m_SamplingSchemeQSpace->GetBThresholdSingleShell()>0)
636  m_SamplingSchemeQSpace->CorrectBValues();
637 
638  // correct values if needed
639  if (m_CorrectDWIValues)
640  CorrectDWI();
641 
642  // normalize values if needed
643  if (m_NormalizeDWI)
644  NormalizeDWI();
645 }
646 
647 template <class TPixelType, unsigned int VImageDimension>
648 void
649 DWIReader<TPixelType, VImageDimension>
650 ::PrintSelf(std::ostream& os, Indent indent) const
651 {
652  Superclass::PrintSelf(os, indent);
653  PrintVar3(true, m_IsInput4DImage, m_NormalizeDWI, m_CorrectDWIValues, os<<indent );
654  PrintVar1(true, m_ConfigurationFile, os<<indent);
655  os << indent << "m_SamplingSchemeQSpace = " << m_SamplingSchemeQSpace << std::endl << std::flush;
656  if (!IsImageEmpty(m_MaskImage))
657  os << indent << "m_MaskImage = " << m_MaskImage << std::endl << std::flush;
658  if (!IsImageEmpty(m_B0Image))
659  os << indent << "m_B0Image = " << m_B0Image << std::endl << std::flush;
660 }
661 
662 }
663 
664 #endif
665 
666 
PrintVar1
#define PrintVar1(cond, var, os)
Definition: utlCoreMacro.h:447
utl.h
helper functions specifically used in dmritool
utl::IsNumber
bool IsNumber(const std::string &input)
Definition: utlCore.h:726
itk::IsImageEmpty
bool IsImageEmpty(const SmartPointer< ImageType > &image)
Definition: utlITK.h:435
itk::DWIReader::PixelType
DWIImageType::PixelType PixelType
Definition: itkDWIReader.h:86
utl::IsFileExist
bool IsFileExist(const std::string &file)
Definition: utlCore.h:529
utl::ReadLinesFirstlineCheck
void ReadLinesFirstlineCheck(const std::string &filename, std::vector< std::vector< std::string > > &strVec, const char *cc=" ")
Definition: utlCore.h:890
itk::DWIReader::PrintSelf
void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
Definition: itkDWIReader.hxx:650
utl::max
const T & max(const T &a, const T &b)
Return the maximum between a and b.
Definition: utlCore.h:263
itk::PrintVariableLengthVector
void PrintVariableLengthVector(const VariableLengthVector< T >vec, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout)
Definition: utlITK.h:45
itk::DWIReader::STDVectorPointer
utl_shared_ptr< STDVectorType > STDVectorPointer
Definition: itkDWIReader.h:81
spams::abs
T abs(const T x)
template version of the fabs function
utlGlobalException
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
itk::VerifyImageInformation
bool VerifyImageInformation(const SmartPointer< Image1Type > &image1, const SmartPointer< Image2Type > &image2, const bool isMinimalDimension=false)
Definition: utlITK.h:850
itk::DWIReader::ReadFromConfigurationFile
void ReadFromConfigurationFile(const std::string &file)
Definition: itkDWIReader.hxx:74
itk::DWIReader::GenerateData
virtual void GenerateData() ITK_OVERRIDE
Definition: itkDWIReader.hxx:628
itk::DWIReader::DetermineIsInput4DImage
static bool DetermineIsInput4DImage(const std::string &dataStr)
Definition: itkDWIReader.hxx:50
utlPrintVar
#define utlPrintVar(cond,...)
Definition: utlCoreMacro.h:309
utl::ReadVector
void ReadVector(const std::string &vectorStr, std::vector< T > &vec, const char *cc=" ")
Definition: utlCore.h:1159
utl::PrintVector
void PrintVector(const std::vector< T > &vec, const std::string &str="", const char *separate=" ", std::ostream &os=std::cout, bool showStats=true)
Definition: utlCore.h:1002
itk::DWIReader::IndexType
DWIImageType::IndexType IndexType
Definition: itkDWIReader.h:84
utl::IsInVector
bool IsInVector(const VectorType &vec, const int size, const T &num, const double eps=1e-10)
Definition: utlCore.h:1394
PrintVar3
#define PrintVar3(cond, var1, var2, var3, os)
Definition: utlCoreMacro.h:462
utlShowPosition
#define utlShowPosition(cond)
Definition: utlCoreMacro.h:554
utl::GetPath
void GetPath(const std::string &fileNameAbsolute, std::string &path, std::string &file)
Definition: utlCore.h:550
itk::DWIReader::MatrixPointer
utl_shared_ptr< MatrixType > MatrixPointer
Definition: itkDWIReader.h:79
Generated by   doxygen 1.8.11