Tomato
itkReadFileListFilter.hxx
1 //
2 // itkReadFileListFilter.h
3 // TomatoLib
4 //
5 // Created by Konrad Werys on 07/11/19.
6 // Copyright © 2017 Konrad Werys. All rights reserved.
7 //
8 
9 #ifndef TomatoLIB_ITKREADFILELISTFILTER_HXX
10 #define TomatoLIB_ITKREADFILELISTFILTER_HXX
11 
12 #include "CmakeConfigForTomato.h"
13 #include "gdcmTomatoReadTags.h"
14 
15 #ifdef USE_ITK
16 
17 #include <iostream>
18 #include <string>
19 #include "itkGDCMImageIO.h"
20 
21 #include "itkStatisticsImageFilter.h"
22 #include "KWUtil.h"
23 #include "itkReadFileListFilter.h"
24 #include "itkDirectory.h"
25 #include "gdcmBase64.h"
26 #include "gdcmTomatoReadTags.h"
27 
28 namespace itk {
29 
30  template< typename TImage >
31  ReadFileListFilter<TImage>::
32  ReadFileListFilter(){
33  m_Verbose = false;
34  }
35 
36  template< typename TImage >
37  void
38  ReadFileListFilter<TImage>
39  ::SetFileList(const std::vector<std::string> fileList){
40  m_FileList = fileList;
41  }
42 
43  template<typename TImage>
44  void ReadFileListFilter<TImage>::SetDirName(const std::string dirName) {
45 
46  m_DirName = dirName;
47  if (m_DirName.empty()) {
48  return;
49  }
50 
51  typename itk::Directory::Pointer directory = itk::Directory::New();
52  directory->Load(m_DirName.c_str());
53 
54  for (size_t i = 0; i < directory->GetNumberOfFiles(); i++){
55  std::string fileName = directory->GetFile(i);
56  if (fileName == "." || fileName == ".." || fileName == ".DS_Store"){
57  continue;
58  }
59  m_FileList.push_back(m_DirName + KWUtil::PathSeparator() + fileName);
60  }
61  }
62 
63  template< typename TImage >
64  template< typename TYPE >
65  vnl_vector<TYPE>
66  ReadFileListFilter<TImage>
67  ::GetVnlVectorFromStdVector(std::vector<TYPE> stdVector){
68  if (stdVector.size() == 0){
69  return vnl_vector<TYPE>();
70  } else {
71  return vnl_vector<TYPE>(&stdVector.front(), stdVector.size());
72  }
73  }
74 
75  template< typename TImage >
76  vnl_vector<double>
77  ReadFileListFilter<TImage>
78  ::GetInvTimes(){
79  if (m_FileList.empty()) return vnl_vector<double>();
80  vnl_vector<double> invTimesVnl = GetVnlVectorFromStdVector(m_InvTimes);
81  vnl_vector<double> invTimesVnl20051572 = GetVnlVectorFromStdVector(m_InvTimes20051572);
82  vnl_vector<double> invTimesVnl00211189 = GetVnlVectorFromStdVector(m_InvTimes00211189);
83  vnl_vector<double> invTimesVnlFromImageCommentsVnl = GetVnlVectorFromStdVector(m_InvTimesFromImageComments);
84  vnl_vector<double> triggerTimesVnl = GetVnlVectorFromStdVector(m_TriggerTimes);
85  if (invTimesVnl.min_value() != invTimesVnl.max_value()) {
86  return invTimesVnl;
87  } else if (invTimesVnl20051572.min_value() != invTimesVnl20051572.max_value()) {
88  return invTimesVnl20051572;
89  } else if (invTimesVnl00211189.min_value() != invTimesVnl00211189.max_value()) {
90  return invTimesVnl00211189;
91  } else if (invTimesVnlFromImageCommentsVnl.min_value() != invTimesVnlFromImageCommentsVnl.max_value()) {
92  return invTimesVnlFromImageCommentsVnl;
93  } else if (triggerTimesVnl.min_value() != triggerTimesVnl.max_value()) {
94  return triggerTimesVnl;
95  } else {
96  return invTimesVnl;
97  }
98  }
99 
100  template< typename TImage >
101  vnl_vector<double>
102  ReadFileListFilter<TImage>
103  ::GetRepTimes(){
104  return GetVnlVectorFromStdVector(m_RepTimes);
105  }
106 
107  template< typename TImage >
108  vnl_vector<double>
109  ReadFileListFilter<TImage>
110  ::GetEchoTimes(){
111  vnl_vector<double> echoTimesVnl = GetVnlVectorFromStdVector(m_EchoTimes);
112  vnl_vector<double> echoTimes00191016Vnl = GetVnlVectorFromStdVector(m_EchoTimes00191016);
113  vnl_vector<double> echoTimes00209158Vnl = GetVnlVectorFromStdVector(m_EchoTimes00209158);
114  vnl_vector<double> echoTimesFromImageCommentsVnl = GetVnlVectorFromStdVector(m_EchoTimesFromImageComments);
115  if (echoTimesVnl.min_value() != echoTimesVnl.max_value()) {
116  return echoTimesVnl;
117  } else if (echoTimesFromImageCommentsVnl.min_value() != echoTimesFromImageCommentsVnl.max_value()) {
118  return echoTimesFromImageCommentsVnl;
119 // } else if (echoTimes00191016Vnl.min_value() != echoTimes00191016Vnl.max_value()) {
120 // return echoTimes00191016Vnl;
121  } else if (echoTimes00209158Vnl.min_value() != echoTimes00209158Vnl.max_value()) {
122  return echoTimes00209158Vnl;
123  } else if (m_EchoTimes.size() == 3 ){
124  // TODO: fix, this works only in siemens T2 3 sample cases
125  vnl_vector<double> temp(3);
126  temp[0] = 0;
127  temp[1] = 25;
128  temp[2] = 55;
129  return temp;
130  } else{
131  return echoTimesVnl;
132  }
133  }
134 
135  template< typename TImage >
136  vnl_vector<double>
137  ReadFileListFilter<TImage>
138  ::GetTriggerTimes(){
139  return GetVnlVectorFromStdVector(m_TriggerTimes);
140  }
141 
142  template< typename TImage >
143  vnl_vector<double>
144  ReadFileListFilter<TImage>
145  ::GetAcqTimes(){
146  return GetVnlVectorFromStdVector(m_AcqTimes);
147  }
148 
149  template< typename TImage >
150  vnl_vector<double>
151  ReadFileListFilter<TImage>
152  ::GetRelAcqTimes(){
153  vnl_vector<double> temp = GetAcqTimes();
154  double mymin = temp.min_value();
155  for (unsigned int i = 0; i < temp.size(); ++i){
156  temp[i] = temp[i] - mymin;
157  }
158  return temp;
159  }
160 
161  template< class TImage >
162  void
163  ReadFileListFilter< TImage >
164  ::GenerateData() {
165  // see https://itk.org/Wiki/ITK/Examples/ImageProcessing/TileImageFilter_CreateVolume
166  typename ImageType3D::Pointer outputImage;
167 
168  m_DicomIO = GDCMImageIO::New();
169  m_DicomIO->LoadPrivateTagsOn();
170 
171  itk::FixedArray< PixelType, 3 > layout;
172  layout[0] = 1;
173  layout[1] = 1;
174  layout[2] = m_FileList.size();
175 
176  // tiler
177  typename TileImageType::Pointer tiler = TileImageType::New();
178  tiler->SetLayout( layout );
179 
180  typename ReaderType::Pointer reader;
181 
182  // the main loop
183  unsigned int inputImageNumber = 0;
184  for ( unsigned int i = 0; i < m_FileList.size(); ++i ){
185 
186  //reader
187  reader = ReaderType::New();
188  reader->SetImageIO(m_DicomIO);
189  reader->SetFileName(m_FileList[i]);
190  m_DicomIO->SetFileName(m_FileList[i]);
191  m_DicomIO->ReadImageInformation();
192  reader->Update();
193 
194  tiler->SetInput(i, reader->GetOutput());
195 
196  try {
197  tiler->Update();
198  //std::cout << std::setprecision(9) << reader->GetOutput()->GetOrigin() << std::endl;
199  m_MetaDataDictionaryArray.push_back(reader->GetMetaDataDictionary());
200  m_InvTimes.push_back(FindInversionTime(reader));
201  m_InvTimes20051572.push_back(FindInversionTime20051572(reader));
202  m_InvTimes00211189.push_back(FindInversionTime00211189(reader));
203  m_InvTimesFromImageComments.push_back(FindInversionTimeFromImageComments(reader));
204  m_RepTimes.push_back(FindRepetitionTime(reader));
205  m_EchoTimes.push_back(FindEchoTime(reader));
206  m_EchoTimes00191016.push_back(FindEchoTime00191016(reader));
207  m_EchoTimes00209158.push_back(FindEchoTime00209158(reader));
208  m_EchoTimesFromImageComments.push_back(FindEchoTimeFromImageComments(reader));
209  m_TriggerTimes.push_back(FindTriggerTime(reader));
210  m_AcqTimes.push_back(FindAcqTime(reader));
211  inputImageNumber++;
212  } catch( itk::ExceptionObject & err ) {
213  std::cerr << "Unable to read file: " << m_FileList[i] << std::endl;
214  std::cerr << "ExceptionObject caught!" << std::endl << err << std::endl;
215  }
216  }
217 
218  if (inputImageNumber == 0) return;
219 
220  tiler->UpdateLargestPossibleRegion();
221  // the final update
222  outputImage = tiler->GetOutput();
223  //outputImage->SetOrigin(reader->GetOutput()->GetOrigin());
224  //outputImage->SetDirection(reader->GetOutput()->GetDirection());
225 
226  try {
227  outputImage->Update();
228  } catch( itk::ExceptionObject & err ) {
229  std::cerr << "ExceptionObject caught!" << std::endl;
230  std::cerr << err << std::endl;
231  }
232 
233  if (GetVerbose()) {
234  //
235  typename ImageType3D::PixelType rescaleIntercept = m_DicomIO->GetRescaleIntercept();
236  typename ImageType3D::PixelType rescaleSlope = m_DicomIO->GetRescaleSlope();
237 
238  typedef itk::StatisticsImageFilter<ImageType3D> StatisticsImageFilterType;
239  typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New();
240  statisticsImageFilter->SetInput(outputImage);
241  statisticsImageFilter->Update();
242 
243  std::cout << "Mean: " << statisticsImageFilter->GetMean();
244  std::cout << ", Std.: " << statisticsImageFilter->GetSigma();
245  std::cout << ", Min: " << statisticsImageFilter->GetMinimum();
246  std::cout << ", Max: " << statisticsImageFilter->GetMaximum();
247 
248  std::cout << ", RescaleIntercept: " << rescaleIntercept << ", RescaleSlope: " << rescaleSlope << std::endl;
249 
250  vcl_cout << "Inversion Times: " << GetInvTimes() << vcl_endl;
251  vcl_cout << "Repetition Times: " << GetRepTimes() << vcl_endl;
252  vcl_cout << "Echo Times: " << GetEchoTimes() << vcl_endl;
253  vcl_cout << "Trigger Times: " << GetTriggerTimes() << vcl_endl;
254  vcl_cout << "Rel Acq Times: " << GetRelAcqTimes() << vcl_endl;
255  }
256 
257  this->GraftOutput(outputImage);
258  }
259 
260  template< class TImage>
261  double
262  ReadFileListFilter< TImage>
263  ::FindInversionTime(ReaderType* reader) {
264 
265  double invTime = 0;
266 
267  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
268  DictionaryType::ConstIterator end = dictionary.End();
269 
270  std::string entryId = "0018|0082";
271  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
272 
273  if ( tagItr != end ) {
274  MetaDataStringType::ConstPointer entryvalueStr =
275  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
276 
277  if (entryvalueStr) {
278  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
279  invTime = KWUtil::StringToNumber<double>(tagvalue);
280  }
281  }
282  return invTime;
283  };
284 
285  template< class TImage>
286  double
287  ReadFileListFilter< TImage>
288  ::FindInversionTime20051572(ReaderType* reader) {
289 
290  double invTime = 0;
291 
292  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
293  DictionaryType::ConstIterator end = dictionary.End();
294 
295  std::string entryId = "2005|1572";
296  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
297 
298  if ( tagItr != end ) {
299  MetaDataStringType::ConstPointer entryvalueStr =
300  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
301 
302  if (entryvalueStr) {
303 // std::string decoded;
304 // std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
305 //
306 // int dlen = gdcm::Base64::GetDecodeLength(tagvalue.c_str(), tagvalue.size() );
307 // decoded.resize(dlen);
308 // gdcm::Base64::Decode(&decoded[0], decoded.size(), tagvalue.c_str(), tagvalue.size());
309 // float_t * fDecoded = (float_t *)&decoded[0];
310 // invTime = (double)(*fDecoded);
311 
312  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
313  invTime = KWUtil::StringToNumber<double>(tagvalue);
314  }
315  }
316  return invTime;
317  };
318 
319  template< class TImage>
320  double
321  ReadFileListFilter< TImage>
322  ::FindInversionTime00211189(ReaderType *reader) {
323 
324  double invTime = 0;
325 
326  std::vector<std::pair<int, int> > tags;
327  tags.push_back(std::pair<int, int>(0x5200, 0x9230));
328  tags.push_back(std::pair<int, int>(0x0021, 0x11fe));
329  tags.push_back(std::pair<int, int>(0x0021, 0x1189));
330 
331  std::string tagvalue;
332  if (gdcmTomatoReadTags(tags, reader->GetFileName(), tagvalue, false) == 0) { // EXIT_SUCCESS
333  invTime = *(double*)&tagvalue[0];
334  }
335 
336  return invTime;
337  };
338 
339 
340  template< class TImage>
341  double
342  ReadFileListFilter< TImage>
343  ::FindInversionTimeFromImageComments(ReaderType* reader) {
344 
345  int invTime = 0;
346 
347  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
348  DictionaryType::ConstIterator end = dictionary.End();
349 
350  std::string entryId = "0020|4000";
351  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
352 
353  if ( tagItr != end ) {
354  MetaDataStringType::ConstPointer entryvalueStr = dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
355 
356  if (entryvalueStr) {
357  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
358  std::sscanf(tagvalue.c_str(), "TIeff %i ms ", &invTime);
359  }
360  }
361  return (double)invTime;
362  };
363 
364  template< class TImage>
365  double
366  ReadFileListFilter< TImage>
367  ::FindRepetitionTime(ReaderType* reader) {
368 
369  double repTime = 0;
370 
371  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
372  DictionaryType::ConstIterator end = dictionary.End();
373 
374  std::string entryId = "0018|0080";
375  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
376  if ( tagItr != end ) {
377  MetaDataStringType::ConstPointer entryvalueStr =
378  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
379 
380  if (entryvalueStr) {
381  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
382  repTime = KWUtil::StringToNumber<double>(tagvalue);
383  }
384  }
385  return repTime;
386  };
387 
388  template< class TImage>
389  double
390  ReadFileListFilter< TImage>
391  ::FindEchoTime(ReaderType* reader) {
392 
393  double echoTime = 0;
394 
395  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
396  DictionaryType::ConstIterator end = dictionary.End();
397 
398  std::string entryId = "0018|0081";
399  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
400  if ( tagItr != end ) {
401  MetaDataStringType::ConstPointer entryvalueStr =
402  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
403 
404  if (entryvalueStr) {
405  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
406  echoTime = KWUtil::StringToNumber<double>(tagvalue);
407  }
408  }
409  return echoTime;
410  };
411 
412  template< class TImage>
413  double
414  ReadFileListFilter< TImage>
415  ::FindEchoTime00191016(ReaderType* reader) {
416 
417  double echoTime = 0;
418 
419  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
420  DictionaryType::ConstIterator end = dictionary.End();
421 
422  std::string entryId = "0019|1016";
423  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
424  if ( tagItr != end ) {
425  MetaDataStringType::ConstPointer entryvalueStr =
426  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
427 
428  if (entryvalueStr) {
429  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
430  echoTime = KWUtil::StringToNumber<double>(tagvalue);
431  }
432  }
433  return echoTime;
434  };
435 
436  template< class TImage>
437  double
438  ReadFileListFilter< TImage>
439  ::FindEchoTime00209158(ReaderType* reader) {
440 
441  int echoTime = 0;
442 
443  std::vector<std::pair<int, int> > tags;
444  tags.push_back(std::pair<int, int>(0x5200, 0x9230));
445  tags.push_back(std::pair<int, int>(0x0020, 0x9111));
446  tags.push_back(std::pair<int, int>(0x0020, 0x9158));
447 
448  std::string tagvalue;
449  if ( gdcmTomatoReadTags(tags, reader->GetFileName(), tagvalue, false) == 0) { // EXIT_SUCCESS
450  std::sscanf(tagvalue.c_str(), "T2 prep. duration = %i ms", &echoTime);
451  }
452 
453  return (double)echoTime;
454  };
455 
456  template< class TImage>
457  double
458  ReadFileListFilter< TImage>
459  ::FindEchoTimeFromImageComments(ReaderType* reader) {
460 
461  int echoTime = 0;
462 
463  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
464  DictionaryType::ConstIterator end = dictionary.End();
465 
466  std::string entryId = "0020|4000";
467  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
468 
469  if ( tagItr != end ) {
470  MetaDataStringType::ConstPointer entryvalueStr = dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
471 
472  if (entryvalueStr) {
473  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
474  std::sscanf(tagvalue.c_str(), "T2 prep. duration = %i ms ", &echoTime);
475  }
476  }
477  return (double)echoTime;
478  };
479 
480  template< class TImage>
481  double
482  ReadFileListFilter< TImage>
483  ::FindTriggerTime(ReaderType* reader) {
484 
485  double triggerTime = 0;
486 
487  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
488  DictionaryType::ConstIterator end = dictionary.End();
489 
490  std::string entryId = "0018|1060";
491  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
492  if ( tagItr != end ) {
493  MetaDataStringType::ConstPointer entryvalueStr =
494  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
495 
496  if (entryvalueStr) {
497  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
498  triggerTime = KWUtil::StringToNumber<double>(tagvalue);
499  }
500  }
501  return triggerTime;
502  };
503 
504  template< class TImage>
505  double
506  ReadFileListFilter< TImage>
507  ::FindAcqTime(ReaderType* reader) {
508 
509  double acqTime = 0;
510 
511  const DictionaryType & dictionary = reader->GetMetaDataDictionary();
512  DictionaryType::ConstIterator end = dictionary.End();
513 
514  std::string entryId = "0008|0032";
515  DictionaryType::ConstIterator tagItr = dictionary.Find(entryId);
516  if ( tagItr != end ) {
517  MetaDataStringType::ConstPointer entryvalueStr =
518  dynamic_cast<const MetaDataStringType *>(tagItr->second.GetPointer());
519 
520  if (entryvalueStr) {
521  std::string tagvalue = entryvalueStr->GetMetaDataObjectValue();
522  acqTime = KWUtil::dicomTime2Seconds<double>(tagvalue);
523  }
524  }
525  return acqTime;
526  }
527 
528  template<typename TImage>
529  const std::vector<std::string> &ReadFileListFilter<TImage>::GetFileList() const {
530  return m_FileList;
531  }
532 
533  template<typename TImage>
534  const std::string &ReadFileListFilter<TImage>::GetDirName() const {
535  return m_DirName;
536  }
537 
538 }
539 
540 #endif // USE_ITK
541 
542 #endif //TomatoLIB_ITKREADFILELISTFILTER_HXX
itk
Definition: itkImageFileReaderKW.h:31