Tomato
Tomato_exportT2.hxx
Go to the documentation of this file.
1 
7 #ifndef TOMATO_TOMATO_EXPORTT2_HXX
8 #define TOMATO_TOMATO_EXPORTT2_HXX
9 
10 #include "CmakeConfigForTomato.h"
11 #ifdef USE_ITK
12 
13 #include "gdcmUIDGenerator.h"
14 #include "itkGDCMImageIO.h"
15 #include "itkGDCMSeriesFileNames.h"
16 #include "itkMultiplyImageFilter.h"
17 #include "itkImageFileWriter.h"
18 #include "itkFileTools.h"
19 #include "itkAdaptImageFilter.h"
20 #include "itkThresholdImageFilter.h"
21 
22 namespace Ox {
23 
24  template< typename MeasureType >
25  int
26  Tomato<MeasureType>
27  ::exportT2ToDicom(){
28 
29  if (_opts->dir_output_map.empty()){
30  printf("No DICOM export, dir_output_map not given.\n");
31  return EXIT_FAILURE;
32  }
33 
34  // some typedefs
35  typedef itk::Colorbar2DImageFilter< ImageType2D > OxColorbarImageFilterType;
36  typedef itk::AdaptImageFilter < ImageType2D, OutputImageType, CastPixelAccessor< InputPixelType, OutputPixelType > > ImageAdaptorType;
37  typedef itk::ImageFileWriter<OutputImageType> WriterType;
38  typedef itk::MultiplyImageFilter< ImageType2D, ImageType2D, ImageType2D > MultiplyImageFilterType;
39  typedef itk::ThresholdImageFilter<ImageType2D> ThresholdImageFilterType;
40 
41 
42  // OxColorbarImageFilter
43  typename OxColorbarImageFilterType::Pointer OxColorbarFilter = OxColorbarImageFilterType::New();
44  OxColorbarFilter->SetInput(_imageCalculatorItk->GetT2Image());
45  OxColorbarFilter->SetAddColorbar(_opts->use_colorbar);
46 
47  // this one I want to reset
48  itk::EncapsulateMetaData<std::string>( _dictionaryInput, std::string("0028|1052"), "0"); // Rescale intercept
49  itk::EncapsulateMetaData<std::string>( _dictionaryInput, std::string("0028|1053"), "1"); // Rescale slope
50 
51  // T2 map with color
52  DictionaryType dictionaryOutput_T2(_dictionaryInput);
53 
54  // UIDs
55  gdcm::UIDGenerator sopuid;
56  gdcm::UIDGenerator suid;
57  std::string sopInstanceUID_T2 = sopuid.Generate();
58  std::string seriesUID_T2;
59  if (_opts->output_map_series_instance_uid.empty()){
60  seriesUID_T2 = suid.Generate();
61  } else {
62  seriesUID_T2 = _opts->output_map_series_instance_uid;
63  }
64 
65  // seriesNumber
66  std::string seriesNumber;
67  itk::ExposeMetaData<std::string>(_dictionaryInput, "0020|0011", seriesNumber);
68  if (_opts->output_map_series_number == 0) {
69  _opts->output_map_series_number = KWUtil::StringToNumber<int>(seriesNumber) + 10002;
70  }
71  std::string newSeriesNumber_T2 = KWUtil::NumberToString(_opts->output_map_series_number);
72 
73  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0008|0018"), sopInstanceUID_T2);
74  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0002|0003"), sopInstanceUID_T2);
75  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0020|000e"), seriesUID_T2);
76 
77  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0020|0011"), newSeriesNumber_T2); // series number
78 // itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0028|1052"), "0"); // rescale intercept
79 // itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0028|1053"), "10"); // rescale slope. Matlab is interpreting it the wrong way
80 
81  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0020|4000"), "T2*10 map"); // series number
82  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2, std::string("0020|0013"), std::string("1")); // instance number
83 
84  // R2 map
85  DictionaryType dictionaryOutput_R2(_dictionaryInput);
86 
87  std::string sopInstanceUID_R2 = sopuid.Generate();
88  std::string seriesUID_R2;
89  if (_opts->output_fitparams_series_instance_uid.empty()){
90  seriesUID_R2 = suid.Generate();
91  } else {
92  seriesUID_R2 = _opts->output_map_series_instance_uid;
93  }
94 
95  itk::ExposeMetaData<std::string>(_dictionaryInput, "0020|0011", seriesNumber);
96  if (_opts->output_fitparams_series_number == 0) {
97  _opts->output_fitparams_series_number = KWUtil::StringToNumber<int>(seriesNumber) + 10003;
98  }
99  std::string newSeriesNumber_R2 = KWUtil::NumberToString(_opts->output_fitparams_series_number);
100 
101 
102  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0008|0018"), sopInstanceUID_R2);
103  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0002|0003"), sopInstanceUID_R2);
104  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0020|000e"), seriesUID_R2);
105 
106  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0020|0011"), newSeriesNumber_R2); // series number
107 
108  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0020|4000"), std::string("Rsquare*4000 Map")); // image comments
109  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0020|0013"), std::string("3")); // instance number
110  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0028|0106"), std::string("0")); // smallest pixel
111  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0028|0107"), std::string("4096")); // largest pixel
112  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0028|1050"), std::string("3900")); // window center
113  itk::EncapsulateMetaData<std::string>( dictionaryOutput_R2, std::string("0028|1051"), std::string("200")); // window width
114 
115  // A map
116  DictionaryType dictionaryOutput_A(_dictionaryInput);
117 
118  std::string sopInstanceUID_A = sopuid.Generate();
119 
120  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0008|0018"), sopInstanceUID_A);
121  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0002|0003"), sopInstanceUID_A);
122  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0020|000e"), seriesUID_R2);
123 
124  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0020|0011"), newSeriesNumber_R2); // series number
125 
126  // itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0020|4000"), std::string("Signal(TI)=[___AMap___]-B*exp(-TI/T2star)"));
127  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0020|4000"), std::string("A")); // image comments
128  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0020|0013"), std::string("4")); // instance number
129  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0028|0106"), std::string("0")); // smallest pixel
130  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0028|0107"), std::string("4096")); // largest pixel
131  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0028|1050"), std::string("2000")); // window center
132  itk::EncapsulateMetaData<std::string>( dictionaryOutput_A, std::string("0028|1051"), std::string("4000")); // window width
133 
134  // B map
135  DictionaryType dictionaryOutput_B(_dictionaryInput);
136 
137  std::string sopInstanceUID_B = sopuid.Generate();
138 
139  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0008|0018"), sopInstanceUID_B);
140  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0002|0003"), sopInstanceUID_B);
141  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0020|000e"), seriesUID_R2);
142 
143  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0020|0011"), newSeriesNumber_R2); // series number
144 
145  // itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0020|4000"), std::string("Signal(TI)=A-[___BMap___]*exp(-TI/T2star)"));
146  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0020|4000"), std::string("B"));
147  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0020|0013"), std::string("5")); // instance number
148  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0028|0106"), std::string("0")); // smallest pixel
149  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0028|0107"), std::string("4096")); // largest pixel
150  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0028|1050"), std::string("2000")); // window center
151  itk::EncapsulateMetaData<std::string>( dictionaryOutput_B, std::string("0028|1051"), std::string("4000")); // window width
152 
153  // T2star map
154  DictionaryType dictionaryOutput_T2star(_dictionaryInput);
155 
156  std::string sopInstanceUID_T2star = sopuid.Generate();
157 
158  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0008|0018"), sopInstanceUID_T2star);
159  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0002|0003"), sopInstanceUID_T2star);
160  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0020|000e"), seriesUID_R2);
161 
162  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0020|0011"), newSeriesNumber_R2); // series number
163 
164  // itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0020|4000"), std::string("Signal(TI)=A-B*exp(-TI/[___T2starMap___])")); // series number
165  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0020|0013"), std::string("6")); // instance number
166  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0028|0106"), std::string("0")); // smallest pixel
167  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0028|0107"), std::string("4096")); // largest pixel
168  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0028|1050"), std::string("2000")); // window center
169  itk::EncapsulateMetaData<std::string>( dictionaryOutput_T2star, std::string("0028|1051"), std::string("4000")); // window width
170 
171  // _RelSNRFinal map
172  DictionaryType dictionaryOutput_RelSNRFinal(_dictionaryInput);
173 
174  std::string sopInstanceUID_RelSNRFinal = sopuid.Generate();
175 
176  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0008|0018"), sopInstanceUID_RelSNRFinal);
177  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0002|0003"), sopInstanceUID_RelSNRFinal);
178  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0020|000e"), seriesUID_R2);
179 
180  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0020|0011"), newSeriesNumber_R2); // series number
181 
182  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0020|4000"), std::string("_RelSNRFinal")); // series number
183  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0020|0013"), std::string("9")); // instance number
184  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0028|0106"), std::string("0")); // smallest pixel
185  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0028|0107"), std::string("4096")); // largest pixel
186  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0028|1050"), std::string("2000")); // window center
187  itk::EncapsulateMetaData<std::string>( dictionaryOutput_RelSNRFinal, std::string("0028|1051"), std::string("4000")); // window width
188 
189  // get GDCMImageIO ready before export
190  typedef itk::GDCMImageIO ImageIOType;
191  ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
192  gdcmImageIO->KeepOriginalUIDOn();
193 
194  // to cast
195  typename ImageAdaptorType::Pointer adaptor = ImageAdaptorType::New();
196 
197  // get the writer ready before export
198  typename WriterType::Pointer writer = WriterType::New();
199  writer->SetUseInputMetaDataDictionary(false);
200 
201  // maybe this? http://www.cplusplus.com/reference/string/string/find_last_of/
202  if (!_opts->dir_output_map.empty()) {
203 
204  printf("Saving to: %s ", _opts->dir_output_map.c_str());
205  //std::cout << "Saving to: " << _opts->dir_output_map << std::flush;
206 
207  itk::FileTools::CreateDirectory(_opts->dir_output_map);
208 
209  // scaling T2 * 10
210  typename MultiplyImageFilterType::Pointer multiplyT2Filter = MultiplyImageFilterType::New();
211  multiplyT2Filter->SetInput(OxColorbarFilter->GetOutput() );
212  multiplyT2Filter->SetConstant( 10. );
213 
214  typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
215  thresholdFilter->SetInput(multiplyT2Filter->GetOutput());
216  // thresholdFilter->SetUpper(pow(2, 16) - 1);
217  thresholdFilter->SetUpper(powl(2, 12) - 1);
218  thresholdFilter->SetOutsideValue(thresholdFilter->GetUpper());
219 
220  // export T2 color
221  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_T2);
222 
223  // export to dicom
224  adaptor->SetInput(thresholdFilter->GetOutput());
225  writer->SetFileName(_opts->dir_output_map + KWUtil::PathSeparator() + newSeriesNumber_T2 + "_T2.dcm");
226  writer->SetInput(adaptor->GetOutput());
227  writer->SetImageIO(gdcmImageIO);
228  try {
229  writer->Update();
230  } catch (itk::ExceptionObject &e) {
231  std::cerr << "Exception in file writer " << std::endl;
232  std::cerr << e << std::endl;
233  }
234 
235  std::cout << " Saved!" << std::endl;
236  }
237 
238  if (!_opts->dir_output_fitparams.empty()) {
239 
240  printf("Saving to: %s ", _opts->dir_output_fitparams.c_str());
241  //std::cout << "Saving to: " << _opts->dir_output_fitparams << std::flush;
242 
243  // scaling R2 * 4000
244  typename MultiplyImageFilterType::Pointer multiplyR2Filter = MultiplyImageFilterType::New();
245  multiplyR2Filter->SetInput(_imageCalculatorItk->GetR2Image() );
246  multiplyR2Filter->SetConstant(4000 );
247 
248  OxColorbarFilter->SetInput(multiplyR2Filter->GetOutput());
249  OxColorbarFilter->SetZerosInsteadOfColorbar(true);
250 
251  // mkdir
252  itk::FileTools::CreateDirectory(_opts->dir_output_fitparams);
253 
254  // export R2
255  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_R2);
256 
257  // export to dicom
258  adaptor->SetInput(OxColorbarFilter->GetOutput());
259  writer->SetFileName(_opts->dir_output_fitparams + KWUtil::PathSeparator() + newSeriesNumber_R2 + "_R2.dcm");
260  writer->SetInput(adaptor->GetOutput());
261  writer->SetImageIO(gdcmImageIO);
262  try {
263  writer->Update();
264  } catch (itk::ExceptionObject &e) {
265  std::cerr << "Exception in file writer " << std::endl;
266  std::cerr << e << std::endl;
267  }
268 
269  // export A
270  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_A);
271 
272  OxColorbarFilter->SetInput(_imageCalculatorItk->GetAImage());
273  OxColorbarFilter->SetZerosInsteadOfColorbar(true);
274 
275  // export to dicom
276  adaptor->SetInput(OxColorbarFilter->GetOutput());
277  writer->SetFileName(_opts->dir_output_fitparams + KWUtil::PathSeparator() + newSeriesNumber_R2 + "_A.dcm");
278  writer->SetInput(adaptor->GetOutput());
279  writer->SetImageIO(gdcmImageIO);
280  try {
281  writer->Update();
282  } catch (itk::ExceptionObject &e) {
283  std::cerr << "Exception in file writer " << std::endl;
284  std::cerr << e << std::endl;
285  }
286 
287  // export B
288  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_B);
289 
290  OxColorbarFilter->SetInput(_imageCalculatorItk->GetBImage());
291  OxColorbarFilter->SetZerosInsteadOfColorbar(true);
292 
293  // export to dicom
294  adaptor->SetInput(OxColorbarFilter->GetOutput());
295  writer->SetFileName(_opts->dir_output_fitparams + KWUtil::PathSeparator() + newSeriesNumber_R2 + "_B.dcm");
296  writer->SetInput(adaptor->GetOutput());
297  writer->SetImageIO(gdcmImageIO);
298  try {
299  writer->Update();
300  } catch (itk::ExceptionObject &e) {
301  std::cerr << "Exception in file writer " << std::endl;
302  std::cerr << e << std::endl;
303  }
304 
305  // export T2star
306  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_T2star);
307 
308  OxColorbarFilter->SetInput(_imageCalculatorItk->GetT1starImage());
309  OxColorbarFilter->SetZerosInsteadOfColorbar(true);
310 
311  // export to dicom
312  adaptor->SetInput(OxColorbarFilter->GetOutput());
313  writer->SetFileName(_opts->dir_output_fitparams + KWUtil::PathSeparator() + newSeriesNumber_R2 + "_T2star.dcm");
314  writer->SetInput(adaptor->GetOutput());
315  writer->SetImageIO(gdcmImageIO);
316  try {
317  writer->Update();
318  } catch (itk::ExceptionObject &e) {
319  std::cerr << "Exception in file writer " << std::endl;
320  std::cerr << e << std::endl;
321  }
322 
323  // TODO: check, values seem a little low
324  // export RelSNRFinal
325  gdcmImageIO->SetMetaDataDictionary(dictionaryOutput_RelSNRFinal);
326 
327  OxColorbarFilter->SetInput(_imageCalculatorItk->GetSNRImage());
328  OxColorbarFilter->SetZerosInsteadOfColorbar(true);
329 
330  // export to dicom
331  adaptor->SetInput(OxColorbarFilter->GetOutput());
332  writer->SetFileName(_opts->dir_output_fitparams + KWUtil::PathSeparator() + newSeriesNumber_R2 + "_RelSNRFinal.dcm");
333  writer->SetInput(adaptor->GetOutput());
334  writer->SetImageIO(gdcmImageIO);
335  try {
336  writer->Update();
337  } catch (itk::ExceptionObject &e) {
338  std::cerr << "Exception in file writer " << std::endl;
339  std::cerr << e << std::endl;
340  }
341  std::cout << " Saved!" << std::endl;
342  }
343 
344  return 0; // EXIT_SUCCESS
345  }
346 
347 } // namespace Ox
348 
349 #endif // USE_ITK
350 
351 #endif //TOMATO_TOMATO_EXPORTT2_HXX
Ox
Definition: OxCalculator.h:19