Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkDicomSeriesReader.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 // uncomment for learning more about the internal sorting mechanisms
14 //#define MBILOG_ENABLE_DEBUG
15 
17 #include <mitkImage.h>
18 #include <mitkImageCast.h>
19 #include <mitkLocaleSwitch.h>
20 
21 #include <itkGDCMSeriesFileNames.h>
22 
23 #include <gdcmAttribute.h>
24 #include <gdcmDirectory.h>
25 #include <gdcmPixmapReader.h>
26 #include <gdcmRAWCodec.h>
27 #include <gdcmScanner.h>
28 #include <gdcmSorter.h>
29 #include <gdcmStringFilter.h>
30 #include <gdcmUIDs.h>
31 
32 #include "mitkProperties.h"
33 
34 namespace mitk
35 {
37  {
38  switch (enumValue)
39  {
41  return "Supported";
43  return "PartlySupported";
45  return "Implemented";
47  return "Unsupported";
48  default:
49  return "<unknown value of enum ReaderImplementationLevel>";
50  };
51  }
52 
54  {
55  switch (enumValue)
56  {
58  return "In Patient";
60  return "At Detector";
62  return "Unknown spacing";
63  default:
64  return "<unknown value of enum PixelSpacingInterpretation>";
65  };
66  }
67 
69  {
70  static bool initialized = false;
71  static TagToPropertyMapType dictionary;
72  if (!initialized)
73  {
74  /*
75  Selection criteria:
76  - no sequences because we cannot represent that
77  - nothing animal related (specied, breed registration number), MITK focusses on human medical image processing.
78  - only general attributes so far
79 
80  When extending this, we should make use of a real dictionary (GDCM/DCMTK and lookup the tag names there)
81  */
82 
83  // Patient module
84  dictionary["0010|0010"] = "dicom.patient.PatientsName";
85  dictionary["0010|0020"] = "dicom.patient.PatientID";
86  dictionary["0010|0030"] = "dicom.patient.PatientsBirthDate";
87  dictionary["0010|0040"] = "dicom.patient.PatientsSex";
88  dictionary["0010|0032"] = "dicom.patient.PatientsBirthTime";
89  dictionary["0010|1000"] = "dicom.patient.OtherPatientIDs";
90  dictionary["0010|1001"] = "dicom.patient.OtherPatientNames";
91  dictionary["0010|2160"] = "dicom.patient.EthnicGroup";
92  dictionary["0010|4000"] = "dicom.patient.PatientComments";
93  dictionary["0012|0062"] = "dicom.patient.PatientIdentityRemoved";
94  dictionary["0012|0063"] = "dicom.patient.DeIdentificationMethod";
95 
96  // General Study module
97  dictionary["0020|000d"] = "dicom.study.StudyInstanceUID";
98  dictionary["0008|0020"] = "dicom.study.StudyDate";
99  dictionary["0008|0030"] = "dicom.study.StudyTime";
100  dictionary["0008|0090"] = "dicom.study.ReferringPhysiciansName";
101  dictionary["0020|0010"] = "dicom.study.StudyID";
102  dictionary["0008|0050"] = "dicom.study.AccessionNumber";
103  dictionary["0008|1030"] = "dicom.study.StudyDescription";
104  dictionary["0008|1048"] = "dicom.study.PhysiciansOfRecord";
105  dictionary["0008|1060"] = "dicom.study.NameOfPhysicianReadingStudy";
106 
107  // General Series module
108  dictionary["0008|0060"] = "dicom.series.Modality";
109  dictionary["0020|000e"] = "dicom.series.SeriesInstanceUID";
110  dictionary["0020|0011"] = "dicom.series.SeriesNumber";
111  dictionary["0020|0060"] = "dicom.series.Laterality";
112  dictionary["0008|0021"] = "dicom.series.SeriesDate";
113  dictionary["0008|0031"] = "dicom.series.SeriesTime";
114  dictionary["0008|1050"] = "dicom.series.PerformingPhysiciansName";
115  dictionary["0018|1030"] = "dicom.series.ProtocolName";
116  dictionary["0008|103e"] = "dicom.series.SeriesDescription";
117  dictionary["0008|1070"] = "dicom.series.OperatorsName";
118  dictionary["0018|0015"] = "dicom.series.BodyPartExamined";
119  dictionary["0018|5100"] = "dicom.series.PatientPosition";
120  dictionary["0028|0108"] = "dicom.series.SmallestPixelValueInSeries";
121  dictionary["0028|0109"] = "dicom.series.LargestPixelValueInSeries";
122 
123  // VOI LUT module
124  dictionary["0028|1050"] = "dicom.voilut.WindowCenter";
125  dictionary["0028|1051"] = "dicom.voilut.WindowWidth";
126  dictionary["0028|1055"] = "dicom.voilut.WindowCenterAndWidthExplanation";
127 
128  // Image Pixel module
129  dictionary["0028|0004"] = "dicom.pixel.PhotometricInterpretation";
130  dictionary["0028|0010"] = "dicom.pixel.Rows";
131  dictionary["0028|0011"] = "dicom.pixel.Columns";
132 
133  // Image Plane module
134  dictionary["0028|0030"] = "dicom.PixelSpacing";
135  dictionary["0018|1164"] = "dicom.ImagerPixelSpacing";
136 
137  // Misc
138  dictionary["0008|0005"] = "dicom.SpecificCharacterSet";
139 
140  initialized = true;
141  }
142 
143  return dictionary;
144  }
145 
147  bool sort,
148  bool check_4d,
149  bool correctTilt,
150  UpdateCallBackMethod callback,
151  Image::Pointer preLoadedImageBlock)
152  {
154 
156  filenames, *node, sort, check_4d, correctTilt, callback, preLoadedImageBlock))
157  {
158  if (filenames.empty())
159  {
160  return nullptr;
161  }
162 
163  return node;
164  }
165  else
166  {
167  return nullptr;
168  }
169  }
170 
172  DataNode &node,
173  bool sort,
174  bool check_4d,
175  bool correctTilt,
176  UpdateCallBackMethod callback,
177  itk::SmartPointer<Image> preLoadedImageBlock)
178  {
179  if (filenames.empty())
180  {
181  MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic.";
182  node.SetData(nullptr);
183  return true; // this is not actually an error but the result is very simple
184  }
185 
186  DcmIoType::Pointer io = DcmIoType::New();
187 
188  try
189  {
190  if (io->CanReadFile(filenames.front().c_str()))
191  {
192  io->SetFileName(filenames.front().c_str());
193  io->ReadImageInformation();
194 
195  if (io->GetPixelType() == itk::ImageIOBase::SCALAR || io->GetPixelType() == itk::ImageIOBase::RGB)
196  {
197  LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock);
198  }
199 
200  if (node.GetData())
201  {
202  return true;
203  }
204  }
205  }
206  catch ( const itk::MemoryAllocationError &e )
207  {
208  MITK_ERROR << "Out of memory. Cannot load DICOM series: " << e.what();
209  }
210  catch ( const std::exception &e )
211  {
212  MITK_ERROR << "Error encountered when loading DICOM series:" << e.what();
213  }
214  catch (...)
215  {
216  MITK_ERROR << "Unspecified error encountered when loading DICOM series.";
217  }
218 
219  return false;
220  }
221 
222  bool DicomSeriesReader::IsDicom(const std::string &filename)
223  {
224  DcmIoType::Pointer io = DcmIoType::New();
225 
226  return io->CanReadFile(filename.c_str());
227  }
228 
229  bool DicomSeriesReader::IsPhilips3DDicom(const std::string &filename)
230  {
231  DcmIoType::Pointer io = DcmIoType::New();
232 
233  if (io->CanReadFile(filename.c_str()))
234  {
235  // Look at header Tag 3001,0010 if it is "Philips3D"
236  gdcm::Reader reader;
237  reader.SetFileName(filename.c_str());
238  reader.Read();
239  gdcm::DataSet &data_set = reader.GetFile().GetDataSet();
240  gdcm::StringFilter sf;
241  sf.SetFile(reader.GetFile());
242 
243  if (data_set.FindDataElement(gdcm::Tag(0x3001, 0x0010)) &&
244  (sf.ToString(gdcm::Tag(0x3001, 0x0010)) == "Philips3D "))
245  {
246  return true;
247  }
248  }
249  return false;
250  }
251 
252  bool DicomSeriesReader::ReadPhilips3DDicom(const std::string &filename, itk::SmartPointer<Image> output_image)
253  {
254  // Now get PhilipsSpecific Tags
255 
256  gdcm::PixmapReader reader;
257  reader.SetFileName(filename.c_str());
258  reader.Read();
259  gdcm::DataSet &data_set = reader.GetFile().GetDataSet();
260  gdcm::StringFilter sf;
261  sf.SetFile(reader.GetFile());
262 
263  gdcm::Attribute<0x0028, 0x0011> dimTagX; // coloumns || sagittal
264  gdcm::Attribute<0x3001, 0x1001, gdcm::VR::UL, gdcm::VM::VM1>
265  dimTagZ; // I have no idea what is VM1. // (Philips specific) // axial
266  gdcm::Attribute<0x0028, 0x0010> dimTagY; // rows || coronal
267  gdcm::Attribute<0x0028, 0x0008> dimTagT; // how many frames
268  gdcm::Attribute<0x0018, 0x602c> spaceTagX; // Spacing in X , unit is "physicalTagx" (usually centimeter)
269  gdcm::Attribute<0x0018, 0x602e> spaceTagY;
270  gdcm::Attribute<0x3001, 0x1003, gdcm::VR::FD, gdcm::VM::VM1> spaceTagZ; // (Philips specific)
271  gdcm::Attribute<0x0018, 0x6024> physicalTagX; // if 3, then spacing params are centimeter
272  gdcm::Attribute<0x0018, 0x6026> physicalTagY;
273  gdcm::Attribute<0x3001, 0x1002, gdcm::VR::US, gdcm::VM::VM1> physicalTagZ; // (Philips specific)
274 
275  dimTagX.Set(data_set);
276  dimTagY.Set(data_set);
277  dimTagZ.Set(data_set);
278  dimTagT.Set(data_set);
279  spaceTagX.Set(data_set);
280  spaceTagY.Set(data_set);
281  spaceTagZ.Set(data_set);
282  physicalTagX.Set(data_set);
283  physicalTagY.Set(data_set);
284  physicalTagZ.Set(data_set);
285 
286  unsigned int dimX = dimTagX.GetValue(), dimY = dimTagY.GetValue(), dimZ = dimTagZ.GetValue(),
287  dimT = dimTagT.GetValue(), physicalX = physicalTagX.GetValue(), physicalY = physicalTagY.GetValue(),
288  physicalZ = physicalTagZ.GetValue();
289 
290  float spaceX = spaceTagX.GetValue(), spaceY = spaceTagY.GetValue(), spaceZ = spaceTagZ.GetValue();
291 
292  if (physicalX == 3) // spacing parameter in cm, have to convert it to mm.
293  spaceX = spaceX * 10;
294 
295  if (physicalY == 3) // spacing parameter in cm, have to convert it to mm.
296  spaceY = spaceY * 10;
297 
298  if (physicalZ == 3) // spacing parameter in cm, have to convert it to mm.
299  spaceZ = spaceZ * 10;
300 
301  // Ok, got all necessary Tags!
302  // Now read Pixeldata (7fe0,0010) X x Y x Z x T Elements
303 
304  const gdcm::Pixmap &pixels = reader.GetPixmap();
305  gdcm::RAWCodec codec;
306 
307  codec.SetPhotometricInterpretation(gdcm::PhotometricInterpretation::MONOCHROME2);
308  codec.SetPixelFormat(pixels.GetPixelFormat());
309  codec.SetPlanarConfiguration(0);
310 
311  gdcm::DataElement out;
312  codec.Decode(data_set.GetDataElement(gdcm::Tag(0x7fe0, 0x0010)), out);
313 
314  const gdcm::ByteValue *bv = out.GetByteValue();
315  const char *new_pixels = bv->GetPointer();
316 
317  // Create MITK Image + Geometry
318  typedef itk::Image<unsigned char, 4>
319  ImageType; // Pixeltype might be different sometimes? Maybe read it out from header
320  ImageType::RegionType myRegion;
321  ImageType::SizeType mySize;
322  ImageType::IndexType myIndex;
323  ImageType::SpacingType mySpacing;
324  ImageType::Pointer imageItk = ImageType::New();
325 
326  mySpacing[0] = spaceX;
327  mySpacing[1] = spaceY;
328  mySpacing[2] = spaceZ;
329  mySpacing[3] = 1;
330  myIndex[0] = 0;
331  myIndex[1] = 0;
332  myIndex[2] = 0;
333  myIndex[3] = 0;
334  mySize[0] = dimX;
335  mySize[1] = dimY;
336  mySize[2] = dimZ;
337  mySize[3] = dimT;
338  myRegion.SetSize(mySize);
339  myRegion.SetIndex(myIndex);
340  imageItk->SetSpacing(mySpacing);
341  imageItk->SetRegions(myRegion);
342  imageItk->Allocate();
343  imageItk->FillBuffer(0);
344 
345  itk::ImageRegionIterator<ImageType> iterator(imageItk, imageItk->GetLargestPossibleRegion());
346  iterator.GoToBegin();
347  unsigned long pixCount = 0;
348  unsigned long planeSize = dimX * dimY;
349  unsigned long planeCount = 0;
350  unsigned long timeCount = 0;
351  unsigned long numberOfSlices = dimZ;
352 
353  while (!iterator.IsAtEnd())
354  {
355  unsigned long adressedPixel =
356  pixCount +
357  (numberOfSlices - 1 - planeCount) * planeSize // add offset to adress the first pixel of current plane
358  + timeCount * numberOfSlices * planeSize; // add time offset
359 
360  iterator.Set(new_pixels[adressedPixel]);
361  pixCount++;
362  ++iterator;
363 
364  if (pixCount == planeSize)
365  {
366  pixCount = 0;
367  planeCount++;
368  }
369  if (planeCount == numberOfSlices)
370  {
371  planeCount = 0;
372  timeCount++;
373  }
374  if (timeCount == dimT)
375  {
376  break;
377  }
378  }
379  mitk::CastToMitkImage(imageItk, output_image);
380  return true; // actually never returns false yet.. but exception possible
381  }
382 
383  std::string DicomSeriesReader::ConstCharStarToString(const char *s) { return s ? std::string(s) : std::string(); }
384  bool DicomSeriesReader::DICOMStringToSpacing(const std::string &s, ScalarType &spacingX, ScalarType &spacingY)
385  {
386  bool successful = false;
387 
388  std::istringstream spacingReader(s);
389  std::string spacing;
390  if (std::getline(spacingReader, spacing, '\\'))
391  {
392  spacingY = atof(spacing.c_str());
393 
394  if (std::getline(spacingReader, spacing, '\\'))
395  {
396  spacingX = atof(spacing.c_str());
397 
398  successful = true;
399  }
400  }
401 
402  return successful;
403  }
404 
405  Point3D DicomSeriesReader::DICOMStringToPoint3D(const std::string &s, bool &successful)
406  {
407  Point3D p;
408  successful = true;
409 
410  std::istringstream originReader(s);
411  std::string coordinate;
412  unsigned int dim(0);
413  while (std::getline(originReader, coordinate, '\\') && dim < 3)
414  {
415  p[dim++] = atof(coordinate.c_str());
416  }
417 
418  if (dim && dim != 3)
419  {
420  successful = false;
421  MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim
422  << " instead of 3 values.";
423  }
424  else if (dim == 0)
425  {
426  successful = false;
427  p.Fill(0.0); // assume default (0,0,0)
428  }
429 
430  return p;
431  }
432 
434  Vector3D &right,
435  Vector3D &up,
436  bool &successful)
437  {
438  successful = true;
439 
440  std::istringstream orientationReader(s);
441  std::string coordinate;
442  unsigned int dim(0);
443  while (std::getline(orientationReader, coordinate, '\\') && dim < 6)
444  {
445  if (dim < 3)
446  {
447  right[dim++] = atof(coordinate.c_str());
448  }
449  else
450  {
451  up[dim++ - 3] = atof(coordinate.c_str());
452  }
453  }
454 
455  if (dim && dim != 6)
456  {
457  successful = false;
458  MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim
459  << " instead of 6 values.";
460  }
461  else if (dim == 0)
462  {
463  // fill with defaults
464  right.Fill(0.0);
465  right[0] = 1.0;
466 
467  up.Fill(0.0);
468  up[1] = 1.0;
469 
470  successful = false;
471  }
472  }
473 
475  const StringContainer &files, bool groupImagesWithGantryTilt, const gdcm::Scanner::MappingType &tagValueMappings_)
476  {
477  // result.first = files that fit ITK's assumption
478  // result.second = files that do not fit, should be run through
479  // AnalyzeFileForITKImageSeriesReaderSpacingAssumption()
480  // again
482 
483  // we const_cast here, because I could not use a map.at(), which would make the code much more readable
484  auto &tagValueMappings = const_cast<gdcm::Scanner::MappingType &>(tagValueMappings_);
485  const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image Position (Patient)
486  const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation
487  const gdcm::Tag tagGantryTilt(0x0018, 0x1120); // gantry tilt
488 
489  Vector3D fromFirstToSecondOrigin;
490  fromFirstToSecondOrigin.Fill(0.0);
491  bool fromFirstToSecondOriginInitialized(false);
492  Point3D thisOrigin;
493  thisOrigin.Fill(0.0f);
494  Point3D lastOrigin;
495  lastOrigin.Fill(0.0f);
496  Point3D lastDifferentOrigin;
497  lastDifferentOrigin.Fill(0.0f);
498 
499  bool lastOriginInitialized(false);
500 
501  MITK_DEBUG << "--------------------------------------------------------------------------------";
502  MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: "
503  << groupImagesWithGantryTilt << ")";
504  unsigned int fileIndex(0);
505  for (auto fileIter = files.begin(); fileIter != files.end(); ++fileIter, ++fileIndex)
506  {
507  bool fileFitsIntoPattern(false);
508  std::string thisOriginString;
509  // Read tag value into point3D. PLEASE replace this by appropriate GDCM code if you figure out how to do that
510  thisOriginString = ConstCharStarToString(tagValueMappings[fileIter->c_str()][tagImagePositionPatient]);
511 
512  if (thisOriginString.empty())
513  {
514  // don't let such files be in a common group. Everything without position information will be loaded as a single
515  // slice:
516  // with standard DICOM files this can happen to: CR, DX, SC
517  MITK_DEBUG << " ==> Sort away " << *fileIter
518  << " for later analysis (no position information)"; // we already have one occupying this position
519 
520  if (result.GetBlockFilenames().empty()) // nothing WITH position information yet
521  {
522  // ==> this is a group of its own, stop processing, come back later
523  result.AddFileToSortedBlock(*fileIter);
524 
525  StringContainer remainingFiles;
526  remainingFiles.insert(remainingFiles.end(), fileIter + 1, files.end());
527  result.AddFilesToUnsortedBlock(remainingFiles);
528 
529  fileFitsIntoPattern = false;
530  break; // no files anymore
531  }
532  else
533  {
534  // ==> this does not match, consider later
535  result.AddFileToUnsortedBlock(*fileIter);
536  fileFitsIntoPattern = false;
537  continue; // next file
538  }
539  }
540 
541  bool ignoredConversionError(-42); // hard to get here, no graceful way to react
542  thisOrigin = DICOMStringToPoint3D(thisOriginString, ignoredConversionError);
543 
544  MITK_DEBUG << " " << fileIndex << " " << *fileIter << " at "
545  /* << thisOriginString */
546  << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")";
547 
548  if (lastOriginInitialized && (thisOrigin == lastOrigin))
549  {
550  MITK_DEBUG << " ==> Sort away " << *fileIter
551  << " for separate time step"; // we already have one occupying this position
552  result.AddFileToUnsortedBlock(*fileIter);
553  fileFitsIntoPattern = false;
554  }
555  else
556  {
557  if (!fromFirstToSecondOriginInitialized &&
558  lastOriginInitialized) // calculate vector as soon as possible when we get a new position
559  {
560  fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin;
561  fromFirstToSecondOriginInitialized = true;
562 
563  // Here we calculate if this slice and the previous one are well aligned,
564  // i.e. we test if the previous origin is on a line through the current
565  // origin, directed into the normal direction of the current slice.
566 
567  // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry,
568  // which cannot be simply loaded into a single mitk::Image at the moment.
569  // For this case, we flag this finding in the result and DicomSeriesReader
570  // can correct for that later.
571 
572  Vector3D right;
573  right.Fill(0.0);
574  Vector3D up;
575  right.Fill(0.0); // might be down as well, but it is just a name at this point
577  tagValueMappings[fileIter->c_str()][tagImageOrientation], right, up, ignoredConversionError);
578 
579  GantryTiltInformation tiltInfo(lastDifferentOrigin, thisOrigin, right, up, 1);
580 
581  if (tiltInfo.IsSheared()) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt
582  {
583  /* optimistic approach, accepting gantry tilt: save file for later, check all further files */
584 
585  // at this point we have TWO slices analyzed! if they are the only two files, we still split, because there
586  // is
587  // no third to verify our tilting assumption.
588  // later with a third being available, we must check if the initial tilting vector is still valid. if yes,
589  // continue.
590  // if NO, we need to split the already sorted part (result.first) and the currently analyzed file
591  // (*fileIter)
592 
593  // tell apart gantry tilt from overall skewedness
594  // sort out irregularly sheared slices, that IS NOT tilting
595 
596  if (groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt())
597  {
598  // check if this is at least roughly the same angle as recorded in DICOM tags
599  if (tagValueMappings[fileIter->c_str()].find(tagGantryTilt) != tagValueMappings[fileIter->c_str()].end())
600  {
601  // read value, compare to calculated angle
602  std::string tiltStr = ConstCharStarToString(tagValueMappings[fileIter->c_str()][tagGantryTilt]);
603  double angle = atof(tiltStr.c_str());
604 
605  MITK_DEBUG << "Comparing recorded tilt angle " << angle << " against calculated value "
606  << tiltInfo.GetTiltAngleInDegrees();
607  // TODO we probably want the signs correct, too (that depends: this is just a rough check, nothing
608  // serious)
609  // TODO TODO TODO when angle -27 and tiltangle 63, this will never trigger the if-clause... useless
610  // check
611  // in this case! old bug..?!
612  if (fabs(angle) - tiltInfo.GetTiltAngleInDegrees() > 0.25)
613  {
614  result.AddFileToUnsortedBlock(*fileIter); // sort away for further analysis
615  fileFitsIntoPattern = false;
616  }
617  else // tilt angle from header is less than 0.25 degrees different from what we calculated, assume this
618  // is
619  // fine
620  {
621  result.FlagGantryTilt();
622  result.AddFileToSortedBlock(*fileIter); // this file is good for current block
623  fileFitsIntoPattern = true;
624  }
625  }
626  else // we cannot check the calculated tilt angle against the one from the dicom header (so we assume we
627  // are
628  // right)
629  {
630  result.FlagGantryTilt();
631  result.AddFileToSortedBlock(*fileIter); // this file is good for current block
632  fileFitsIntoPattern = true;
633  }
634  }
635  else // caller does not want tilt compensation OR shearing is more complicated than tilt
636  {
637  result.AddFileToUnsortedBlock(*fileIter); // sort away for further analysis
638  fileFitsIntoPattern = false;
639  }
640  }
641  else // not sheared
642  {
643  result.AddFileToSortedBlock(*fileIter); // this file is good for current block
644  fileFitsIntoPattern = true;
645  }
646  }
647  else if (fromFirstToSecondOriginInitialized) // we already know the offset between slices
648  {
649  Point3D assumedOrigin = lastDifferentOrigin + fromFirstToSecondOrigin;
650 
651  Vector3D originError = assumedOrigin - thisOrigin;
652  double norm = originError.GetNorm();
653  double toleratedError(0.005); // max. 1/10mm error when measurement crosses 20 slices in z direction
654 
655  if (norm > toleratedError)
656  {
657  MITK_DEBUG << " File does not fit into the inter-slice distance pattern (diff = " << norm << ", allowed "
658  << toleratedError << ").";
659  MITK_DEBUG << " Expected position (" << assumedOrigin[0] << "," << assumedOrigin[1] << ","
660  << assumedOrigin[2] << "), got position (" << thisOrigin[0] << "," << thisOrigin[1] << ","
661  << thisOrigin[2] << ")";
662  MITK_DEBUG << " ==> Sort away " << *fileIter << " for later analysis";
663 
664  // At this point we know we deviated from the expectation of ITK's ImageSeriesReader
665  // We split the input file list at this point, i.e. all files up to this one (excluding it)
666  // are returned as group 1, the remaining files (including the faulty one) are group 2
667 
668  /* Optimistic approach: check if any of the remaining slices fits in */
669  result.AddFileToUnsortedBlock(*fileIter); // sort away for further analysis
670  fileFitsIntoPattern = false;
671  }
672  else
673  {
674  result.AddFileToSortedBlock(*fileIter); // this file is good for current block
675  fileFitsIntoPattern = true;
676  }
677  }
678  else // this should be the very first slice
679  {
680  result.AddFileToSortedBlock(*fileIter); // this file is good for current block
681  fileFitsIntoPattern = true;
682  }
683  }
684 
685  // record current origin for reference in later iterations
686  if (!lastOriginInitialized || (fileFitsIntoPattern && (thisOrigin != lastOrigin)))
687  {
688  lastDifferentOrigin = thisOrigin;
689  }
690 
691  lastOrigin = thisOrigin;
692  lastOriginInitialized = true;
693  }
694 
695  if (result.ContainsGantryTilt())
696  {
697  // check here how many files were grouped.
698  // IF it was only two files AND we assume tiltedness (e.g. save "distance")
699  // THEN we would want to also split the two previous files (simple) because
700  // we don't have any reason to assume they belong together
701 
702  if (result.GetBlockFilenames().size() == 2)
703  {
704  result.UndoPrematureGrouping();
705  }
706  }
707 
708  return result;
709  }
710 
712  bool groupImagesWithGantryTilt,
713  const StringContainer &restrictions)
714  {
715  return GetSeries(files, true, groupImagesWithGantryTilt, restrictions);
716  }
717 
719  bool sortTo3DPlust,
720  bool groupImagesWithGantryTilt,
721  const StringContainer & /*restrictions*/)
722  {
733  // use GDCM directly, itk::GDCMSeriesFileNames does not work with GDCM 2
734 
735  // PART I: scan files for sorting relevant DICOM tags,
736  // separate images that differ in any of those
737  // attributes (they cannot possibly form a 3D block)
738 
739  // scan for relevant tags in dicom files
740  gdcm::Scanner scanner;
741  const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP class UID
742  scanner.AddTag(tagSOPClassUID);
743 
744  const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID
745  scanner.AddTag(tagSeriesInstanceUID);
746 
747  const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation
748  scanner.AddTag(tagImageOrientation);
749 
750  const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing
751  scanner.AddTag(tagPixelSpacing);
752 
753  const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing
754  scanner.AddTag(tagImagerPixelSpacing);
755 
756  const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness
757  scanner.AddTag(tagSliceThickness);
758 
759  const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows
760  scanner.AddTag(tagNumberOfRows);
761 
762  const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols
763  scanner.AddTag(tagNumberOfColumns);
764 
765  const gdcm::Tag tagGantryTilt(0x0018, 0x1120); // gantry tilt
766  scanner.AddTag(tagGantryTilt);
767 
768  const gdcm::Tag tagModality(0x0008, 0x0060); // modality
769  scanner.AddTag(tagModality);
770 
771  const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames
772  scanner.AddTag(tagNumberOfFrames);
773 
774  // additional tags read in this scan to allow later analysis
775  // THESE tag are not used for initial separating of files
776  const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image Position (Patient)
777  scanner.AddTag(tagImagePositionPatient);
778 
779  // TODO add further restrictions from arguments (when anybody asks for it)
780 
781  FileNamesGrouping result;
782 
783  // let GDCM scan files
784  if (!scanner.Scan(files))
785  {
786  MITK_ERROR << "gdcm::Scanner failed when scanning " << files.size() << " input files.";
787  return result;
788  }
789 
790  // assign files IDs that will separate them for loading into image blocks
791  for (auto fileIter = scanner.Begin(); fileIter != scanner.End(); ++fileIter)
792  {
793  if (std::string(fileIter->first).empty())
794  continue; // TODO understand why Scanner has empty string entries
795  if (std::string(fileIter->first) == std::string("DICOMDIR"))
796  continue;
797 
798  /* sort out multi-frame
799  if ( scanner.GetValue( fileIter->first , tagNumberOfFrames ) )
800  {
801  MITK_INFO << "Ignoring " << fileIter->first << " because we cannot handle multi-frame images.";
802  continue;
803  }
804  */
805 
806  // we const_cast here, because I could not use a map.at() function in CreateMoreUniqueSeriesIdentifier.
807  // doing the same thing with find would make the code less readable. Since we forget the Scanner results
808  // anyway after this function, we can simply tolerate empty map entries introduced by bad operator[] access
809  std::string moreUniqueSeriesId =
810  CreateMoreUniqueSeriesIdentifier(const_cast<gdcm::Scanner::TagToValue &>(fileIter->second));
811  result[moreUniqueSeriesId].AddFile(fileIter->first);
812  }
813 
814  // PART II: sort slices spatially (or at least consistently if this is NOT possible, see method)
815 
816  for (FileNamesGrouping::const_iterator groupIter = result.begin(); groupIter != result.end(); ++groupIter)
817  {
818  try
819  {
820  result[groupIter->first] =
821  ImageBlockDescriptor(SortSeriesSlices(groupIter->second.GetFilenames())); // sort each slice group spatially
822  }
823  catch (...)
824  {
825  MITK_ERROR << "Caught something.";
826  }
827  }
828 
829  // PART III: analyze pre-sorted images for valid blocks (i.e. blocks of equal z-spacing),
830  // separate into multiple blocks if necessary.
831  //
832  // Analysis performs the following steps:
833  // * imitate itk::ImageSeriesReader: use the distance between the first two images as z-spacing
834  // * check what images actually fulfill ITK's z-spacing assumption
835  // * separate all images that fail the test into new blocks, re-iterate analysis for these blocks
836  // * this includes images which DO NOT PROVIDE spatial information, i.e. all images w/o
837  // ImagePositionPatient will be loaded separately
838 
839  FileNamesGrouping groupsOf3DPlusTBlocks; // final result of this function
840  for (FileNamesGrouping::const_iterator groupIter = result.begin(); groupIter != result.end(); ++groupIter)
841  {
842  FileNamesGrouping groupsOf3DBlocks; // intermediate result for only this group(!)
843  StringContainer filesStillToAnalyze = groupIter->second.GetFilenames();
844  std::string groupUID = groupIter->first;
845  unsigned int subgroup(0);
846  MITK_DEBUG << "Analyze group " << groupUID << " of " << groupIter->second.GetFilenames().size() << " files";
847 
848  while (!filesStillToAnalyze.empty()) // repeat until all files are grouped somehow
849  {
851  filesStillToAnalyze, groupImagesWithGantryTilt, scanner.GetMappings());
852 
853  // enhance the UID for additional groups
854  std::stringstream newGroupUID;
855  newGroupUID << groupUID << '.' << subgroup;
856 
857  ImageBlockDescriptor thisBlock(analysisResult.GetBlockFilenames());
858 
859  std::string firstFileInBlock = thisBlock.GetFilenames().front();
860 
861  thisBlock.SetImageBlockUID(newGroupUID.str());
862  thisBlock.SetSeriesInstanceUID(
863  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagSeriesInstanceUID)));
864  thisBlock.SetHasGantryTiltCorrected(analysisResult.ContainsGantryTilt());
865  thisBlock.SetSOPClassUID(
866  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagSOPClassUID)));
867  thisBlock.SetNumberOfFrames(
868  ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagNumberOfFrames)));
869  thisBlock.SetModality(
870  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagModality)));
871  thisBlock.SetPixelSpacingInformation(
872  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagPixelSpacing)),
873  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(firstFileInBlock.c_str(), tagImagerPixelSpacing)));
874  thisBlock.SetHasMultipleTimePoints(false);
875 
876  groupsOf3DBlocks[newGroupUID.str()] = thisBlock;
877 
878  // MITK_DEBUG << "Result: sorted 3D group " << newGroupUID.str() << " with " << groupsOf3DBlocks[
879  // newGroupUID.str() ].GetFilenames().size() << " files";
880  MITK_DEBUG << "Result: sorted 3D group with " << groupsOf3DBlocks[newGroupUID.str()].GetFilenames().size()
881  << " files";
882  StringContainer debugOutputFiles = analysisResult.GetBlockFilenames();
883  for (StringContainer::const_iterator siter = debugOutputFiles.begin(); siter != debugOutputFiles.end(); ++siter)
884  MITK_DEBUG << " IN " << *siter;
885 
886  ++subgroup;
887 
888  filesStillToAnalyze = analysisResult.GetUnsortedFilenames(); // remember what needs further analysis
889  for (StringContainer::const_iterator siter = filesStillToAnalyze.begin(); siter != filesStillToAnalyze.end();
890  ++siter)
891  MITK_DEBUG << " OUT " << *siter;
892  }
893 
894  // end of grouping, now post-process groups
895 
896  // PART IV: attempt to group blocks to 3D+t blocks if requested
897  // inspect entries of groupsOf3DBlocks
898  // - if number of files is identical to previous entry, collect for 3D+t block
899  // - as soon as number of files changes from previous entry, record collected blocks as 3D+t block,
900  // start
901  // a new one, continue
902 
903  // decide whether or not to group 3D blocks into 3D+t blocks where possible
904  if (!sortTo3DPlust)
905  {
906  // copy 3D blocks to output
907  groupsOf3DPlusTBlocks.insert(groupsOf3DBlocks.begin(), groupsOf3DBlocks.end());
908  }
909  else
910  {
911  // sort 3D+t (as described in "PART IV")
912 
913  MITK_DEBUG << "================================================================================";
914  MITK_DEBUG << "3D+t analysis:";
915 
916  unsigned int numberOfFilesInPreviousBlock(0);
917  std::string previousBlockKey;
918 
919  for (FileNamesGrouping::const_iterator block3DIter = groupsOf3DBlocks.begin();
920  block3DIter != groupsOf3DBlocks.end();
921  ++block3DIter)
922  {
923  unsigned int numberOfFilesInThisBlock = block3DIter->second.GetFilenames().size();
924  std::string thisBlockKey = block3DIter->first;
925 
926  if (numberOfFilesInPreviousBlock == 0)
927  {
928  numberOfFilesInPreviousBlock = numberOfFilesInThisBlock;
929  groupsOf3DPlusTBlocks[thisBlockKey] = block3DIter->second;
930  MITK_DEBUG << " 3D+t group " << thisBlockKey;
931  previousBlockKey = thisBlockKey;
932  }
933  else
934  {
935  bool identicalOrigins;
936  try
937  {
938  // check whether this and the previous block share a comon origin
939  // TODO should be safe, but a little try/catch or other error handling wouldn't hurt
940 
941  const char *origin_value = scanner.GetValue(groupsOf3DBlocks[thisBlockKey].GetFilenames().front().c_str(),
942  tagImagePositionPatient),
943  *previous_origin_value = scanner.GetValue(
944  groupsOf3DBlocks[previousBlockKey].GetFilenames().front().c_str(), tagImagePositionPatient),
945  *destination_value = scanner.GetValue(
946  groupsOf3DBlocks[thisBlockKey].GetFilenames().back().c_str(), tagImagePositionPatient),
947  *previous_destination_value = scanner.GetValue(
948  groupsOf3DBlocks[previousBlockKey].GetFilenames().back().c_str(), tagImagePositionPatient);
949 
950  if (!origin_value || !previous_origin_value || !destination_value || !previous_destination_value)
951  {
952  identicalOrigins = false;
953  }
954  else
955  {
956  std::string thisOriginString = ConstCharStarToString(origin_value);
957  std::string previousOriginString = ConstCharStarToString(previous_origin_value);
958 
959  // also compare last origin, because this might differ if z-spacing is different
960  std::string thisDestinationString = ConstCharStarToString(destination_value);
961  std::string previousDestinationString = ConstCharStarToString(previous_destination_value);
962 
963  identicalOrigins =
964  ((thisOriginString == previousOriginString) && (thisDestinationString == previousDestinationString));
965  }
966  }
967  catch (...)
968  {
969  identicalOrigins = false;
970  }
971 
972  if (identicalOrigins && (numberOfFilesInPreviousBlock == numberOfFilesInThisBlock))
973  {
974  // group with previous block
975  groupsOf3DPlusTBlocks[previousBlockKey].AddFiles(block3DIter->second.GetFilenames());
976  groupsOf3DPlusTBlocks[previousBlockKey].SetHasMultipleTimePoints(true);
977  MITK_DEBUG << " --> group enhanced with another timestep";
978  }
979  else
980  {
981  // start a new block
982  groupsOf3DPlusTBlocks[thisBlockKey] = block3DIter->second;
983  int numberOfTimeSteps =
984  groupsOf3DPlusTBlocks[previousBlockKey].GetFilenames().size() / numberOfFilesInPreviousBlock;
985  MITK_DEBUG << " ==> group closed with " << numberOfTimeSteps << " time steps";
986  previousBlockKey = thisBlockKey;
987  MITK_DEBUG << " 3D+t group " << thisBlockKey << " started";
988  }
989  }
990 
991  numberOfFilesInPreviousBlock = numberOfFilesInThisBlock;
992  }
993  }
994  }
995 
996  MITK_DEBUG << "================================================================================";
997  MITK_DEBUG << "Summary: ";
998  for (FileNamesGrouping::const_iterator groupIter = groupsOf3DPlusTBlocks.begin();
999  groupIter != groupsOf3DPlusTBlocks.end();
1000  ++groupIter)
1001  {
1002  ImageBlockDescriptor block = groupIter->second;
1003  MITK_DEBUG << " " << block.GetFilenames().size() << " '" << block.GetModality() << "' images ("
1004  << block.GetSOPClassUIDAsString() << ") in volume " << block.GetImageBlockUID();
1005  MITK_DEBUG << " (gantry tilt : " << (block.HasGantryTiltCorrected() ? "Yes" : "No") << "; "
1006  "pixel spacing : "
1008  "3D+t: "
1009  << (block.HasMultipleTimePoints() ? "Yes" : "No") << "; "
1010  "reader support: "
1012  StringContainer debugOutputFiles = block.GetFilenames();
1013  for (StringContainer::const_iterator siter = debugOutputFiles.begin(); siter != debugOutputFiles.end(); ++siter)
1014  MITK_DEBUG << " F " << *siter;
1015  }
1016  MITK_DEBUG << "================================================================================";
1017 
1018  return groupsOf3DPlusTBlocks;
1019  }
1020 
1022  bool groupImagesWithGantryTilt,
1023  const StringContainer &restrictions)
1024  {
1025  gdcm::Directory directoryLister;
1026  directoryLister.Load(dir.c_str(), false); // non-recursive
1027  return GetSeries(directoryLister.GetFilenames(), groupImagesWithGantryTilt, restrictions);
1028  }
1029 
1030  std::string DicomSeriesReader::CreateSeriesIdentifierPart(gdcm::Scanner::TagToValue &tagValueMap,
1031  const gdcm::Tag &tag)
1032  {
1033  std::string result;
1034  try
1035  {
1036  result = IDifyTagValue(tagValueMap[tag] ? tagValueMap[tag] : std::string(""));
1037  }
1038  catch ( const std::exception & )
1039  {
1040  // we are happy with even nothing, this will just group images of a series
1041  // MITK_WARN << "Could not access tag " << tag << ": " << e.what();
1042  }
1043 
1044  return result;
1045  }
1046 
1047  std::string DicomSeriesReader::CreateMoreUniqueSeriesIdentifier(gdcm::Scanner::TagToValue &tagValueMap)
1048  {
1049  const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID
1050  const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation
1051  const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing
1052  const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing
1053  const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness
1054  const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows
1055  const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols
1056  const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames
1057 
1058  const char *tagSeriesInstanceUid = tagValueMap[tagSeriesInstanceUID];
1059  if (!tagSeriesInstanceUid)
1060  {
1061  mitkThrow() << "CreateMoreUniqueSeriesIdentifier() could not access series instance UID. Something is seriously "
1062  "wrong with this image, so stopping here.";
1063  }
1064  std::string constructedID = tagSeriesInstanceUid;
1065 
1066  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagNumberOfRows);
1067  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagNumberOfColumns);
1068  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagPixelSpacing);
1069  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagImagerPixelSpacing);
1070  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagSliceThickness);
1071  constructedID += CreateSeriesIdentifierPart(tagValueMap, tagNumberOfFrames);
1072 
1073  // be a bit tolerant for orienatation, let only the first few digits matter
1074  // (http://bugs.mitk.org/show_bug.cgi?id=12263)
1075  // NOT constructedID += CreateSeriesIdentifierPart( tagValueMap, tagImageOrientation );
1076  if (tagValueMap.find(tagImageOrientation) != tagValueMap.end())
1077  {
1078  bool conversionError(false);
1079  Vector3D right;
1080  right.Fill(0.0);
1081  Vector3D up;
1082  right.Fill(0.0);
1083  DICOMStringToOrientationVectors(tagValueMap[tagImageOrientation], right, up, conversionError);
1084  // string newstring sprintf(simplifiedOrientationString, "%.3f\\%.3f\\%.3f\\%.3f\\%.3f\\%.3f", right[0],
1085  // right[1],
1086  // right[2], up[0], up[1], up[2]);
1087  std::ostringstream ss;
1088  ss.setf(std::ios::fixed, std::ios::floatfield);
1089  ss.precision(5);
1090  ss << right[0] << "\\" << right[1] << "\\" << right[2] << "\\" << up[0] << "\\" << up[1] << "\\" << up[2];
1091  std::string simplifiedOrientationString(ss.str());
1092 
1093  constructedID += IDifyTagValue(simplifiedOrientationString);
1094  }
1095 
1096  constructedID.resize(constructedID.length() - 1); // cut of trailing '.'
1097 
1098  return constructedID;
1099  }
1100 
1101  std::string DicomSeriesReader::IDifyTagValue(const std::string &value)
1102  {
1103  std::string IDifiedValue(value);
1104  if (value.empty())
1105  throw std::logic_error("IDifyTagValue() illegaly called with empty tag value");
1106 
1107  // Eliminate non-alnum characters, including whitespace...
1108  // that may have been introduced by concats.
1109  for (std::size_t i = 0; i < IDifiedValue.size(); i++)
1110  {
1111  while (i < IDifiedValue.size() &&
1112  !(IDifiedValue[i] == '.' || (IDifiedValue[i] >= 'a' && IDifiedValue[i] <= 'z') ||
1113  (IDifiedValue[i] >= '0' && IDifiedValue[i] <= '9') ||
1114  (IDifiedValue[i] >= 'A' && IDifiedValue[i] <= 'Z')))
1115  {
1116  IDifiedValue.erase(i, 1);
1117  }
1118  }
1119 
1120  IDifiedValue += ".";
1121  return IDifiedValue;
1122  }
1123 
1125  const std::string &series_uid,
1126  bool groupImagesWithGantryTilt,
1127  const StringContainer &restrictions)
1128  {
1129  FileNamesGrouping allSeries = GetSeries(dir, groupImagesWithGantryTilt, restrictions);
1130  StringContainer resultingFileList;
1131 
1132  for (FileNamesGrouping::const_iterator idIter = allSeries.begin(); idIter != allSeries.end(); ++idIter)
1133  {
1134  if (idIter->first.find(series_uid) == 0) // this ID starts with given series_uid
1135  {
1136  return idIter->second.GetFilenames();
1137  }
1138  }
1139 
1140  return resultingFileList;
1141  }
1142 
1144  {
1145  /* we CAN expect a group of equal
1146  - series instance uid
1147  - image orientation
1148  - pixel spacing
1149  - imager pixel spacing
1150  - slice thickness
1151  - number of rows/columns
1152 
1153  (each piece of information except the rows/columns might be missing)
1154 
1155  sorting with GdcmSortFunction tries its best by sorting by spatial position
1156  and more hints (acquisition number, acquisition time, trigger time) but will
1157  always produce a sorting by falling back to SOP Instance UID.
1158  */
1159  gdcm::Sorter sorter;
1160 
1161  sorter.SetSortFunction(DicomSeriesReader::GdcmSortFunction);
1162  try
1163  {
1164  if (sorter.Sort(unsortedFilenames))
1165  {
1166  return sorter.GetFilenames();
1167  }
1168  else
1169  {
1170  MITK_WARN << "Sorting error. Leaving series unsorted.";
1171  return unsortedFilenames;
1172  }
1173  }
1174  catch ( const std::logic_error & )
1175  {
1176  MITK_WARN << "Sorting error. Leaving series unsorted.";
1177  return unsortedFilenames;
1178  }
1179  }
1180 
1181  bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2)
1182  {
1183  // This method MUST accept missing location and position information (and all else, too)
1184  // because we cannot rely on anything
1185  // (restriction on the sentence before: we have to provide consistent sorting, so we
1186  // rely on the minimum information all DICOM files need to provide: SOP Instance UID)
1187 
1188  /* we CAN expect a group of equal
1189  - series instance uid
1190  - image orientation
1191  - pixel spacing
1192  - imager pixel spacing
1193  - slice thickness
1194  - number of rows/columns
1195  */
1196  static const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image Position (Patient)
1197  static const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation
1198 
1199  // see if we have Image Position and Orientation
1200  if (ds1.FindDataElement(tagImagePositionPatient) && ds1.FindDataElement(tagImageOrientation) &&
1201  ds2.FindDataElement(tagImagePositionPatient) && ds2.FindDataElement(tagImageOrientation))
1202  {
1203  gdcm::Attribute<0x0020, 0x0032> image_pos1; // Image Position (Patient)
1204  gdcm::Attribute<0x0020, 0x0037> image_orientation1; // Image Orientation (Patient)
1205 
1206  image_pos1.Set(ds1);
1207  image_orientation1.Set(ds1);
1208 
1209  gdcm::Attribute<0x0020, 0x0032> image_pos2;
1210  gdcm::Attribute<0x0020, 0x0037> image_orientation2;
1211 
1212  image_pos2.Set(ds2);
1213  image_orientation2.Set(ds2);
1214 
1215  /*
1216  we tolerate very small differences in image orientation, since we got to know about
1217  acquisitions where these values change across a single series (7th decimal digit)
1218  (http://bugs.mitk.org/show_bug.cgi?id=12263)
1219 
1220  still, we want to check if our assumption of 'almost equal' orientations is valid
1221  */
1222  for (unsigned int dim = 0; dim < 6; ++dim)
1223  {
1224  if (fabs(image_orientation2[dim] - image_orientation1[dim]) > 0.0001)
1225  {
1226  MITK_ERROR << "Dicom images have different orientations.";
1227  throw std::logic_error(
1228  "Dicom images have different orientations. Call GetSeries() first to separate images.");
1229  }
1230  }
1231 
1232  double normal[3];
1233 
1234  normal[0] = image_orientation1[1] * image_orientation1[5] - image_orientation1[2] * image_orientation1[4];
1235  normal[1] = image_orientation1[2] * image_orientation1[3] - image_orientation1[0] * image_orientation1[5];
1236  normal[2] = image_orientation1[0] * image_orientation1[4] - image_orientation1[1] * image_orientation1[3];
1237 
1238  double dist1 = 0.0, dist2 = 0.0;
1239 
1240  // this computes the distance from world origin (0,0,0) ALONG THE NORMAL of the image planes
1241  for (unsigned char i = 0u; i < 3u; ++i)
1242  {
1243  dist1 += normal[i] * image_pos1[i];
1244  dist2 += normal[i] * image_pos2[i];
1245  }
1246 
1247  // if we can sort by just comparing the distance, we do exactly that
1248  if (fabs(dist1 - dist2) >= mitk::eps)
1249  {
1250  // default: compare position
1251  return dist1 < dist2;
1252  }
1253  else // we need to check more properties to distinguish slices
1254  {
1255  // try to sort by Acquisition Number
1256  static const gdcm::Tag tagAcquisitionNumber(0x0020, 0x0012);
1257  if (ds1.FindDataElement(tagAcquisitionNumber) && ds2.FindDataElement(tagAcquisitionNumber))
1258  {
1259  gdcm::Attribute<0x0020, 0x0012> acquisition_number1; // Acquisition number
1260  gdcm::Attribute<0x0020, 0x0012> acquisition_number2;
1261 
1262  acquisition_number1.Set(ds1);
1263  acquisition_number2.Set(ds2);
1264 
1265  if (acquisition_number1 != acquisition_number2)
1266  {
1267  return acquisition_number1 < acquisition_number2;
1268  }
1269  else // neither position nor acquisition number are good for sorting, so check more
1270  {
1271  // try to sort by Acquisition Time
1272  static const gdcm::Tag tagAcquisitionTime(0x0008, 0x0032);
1273  if (ds1.FindDataElement(tagAcquisitionTime) && ds2.FindDataElement(tagAcquisitionTime))
1274  {
1275  gdcm::Attribute<0x0008, 0x0032> acquisition_time1; // Acquisition time
1276  gdcm::Attribute<0x0008, 0x0032> acquisition_time2;
1277 
1278  acquisition_time1.Set(ds1);
1279  acquisition_time2.Set(ds2);
1280 
1281  if (acquisition_time1 != acquisition_time2)
1282  {
1283  return acquisition_time1 < acquisition_time2;
1284  }
1285  else // we gave up on image position, acquisition number and acquisition time now
1286  {
1287  // let's try trigger time
1288  static const gdcm::Tag tagTriggerTime(0x0018, 0x1060);
1289  if (ds1.FindDataElement(tagTriggerTime) && ds2.FindDataElement(tagTriggerTime))
1290  {
1291  gdcm::Attribute<0x0018, 0x1060> trigger_time1; // Trigger time
1292  gdcm::Attribute<0x0018, 0x1060> trigger_time2;
1293 
1294  trigger_time1.Set(ds1);
1295  trigger_time2.Set(ds2);
1296 
1297  if (trigger_time1 != trigger_time2)
1298  {
1299  return trigger_time1 < trigger_time2;
1300  }
1301  // ELSE!
1302  // for this and many previous ifs we fall through if nothing lets us sort
1303  } // .
1304  } // .
1305  } // .
1306  }
1307  }
1308  }
1309  } // .
1310 
1311  // LAST RESORT: all valuable information for sorting is missing.
1312  // Sort by some meaningless but unique identifiers to satisfy the sort function
1313  static const gdcm::Tag tagSOPInstanceUID(0x0008, 0x0018);
1314  if (ds1.FindDataElement(tagSOPInstanceUID) && ds2.FindDataElement(tagSOPInstanceUID))
1315  {
1316  MITK_DEBUG
1317  << "Dicom images are missing attributes for a meaningful sorting, falling back to SOP instance UID comparison.";
1318  gdcm::Attribute<0x0008, 0x0018> SOPInstanceUID1; // SOP instance UID is mandatory and unique
1319  gdcm::Attribute<0x0008, 0x0018> SOPInstanceUID2;
1320 
1321  SOPInstanceUID1.Set(ds1);
1322  SOPInstanceUID2.Set(ds2);
1323 
1324  return SOPInstanceUID1 < SOPInstanceUID2;
1325  }
1326  else
1327  {
1328  // no DICOM file should really come down here, this should only be reached with unskillful and unlucky
1329  // manipulation
1330  // of files
1331  std::string error_message("Malformed DICOM images, which do not even contain a SOP Instance UID.");
1332  MITK_ERROR << error_message;
1333  throw std::logic_error(error_message);
1334  }
1335  }
1336 
1338  {
1339  std::stringstream configuration;
1340  configuration << "MITK_USE_GDCMIO: ";
1341  configuration << "true";
1342  configuration << "\n";
1343 
1344  configuration << "GDCM_VERSION: ";
1345 #ifdef GDCM_MAJOR_VERSION
1346  configuration << GDCM_VERSION;
1347 #endif
1348  // configuration << "\n";
1349 
1350  return configuration.str();
1351  }
1352 
1354  const gdcm::Scanner::MappingType &tagValueMappings_,
1355  DcmIoType *io,
1356  const ImageBlockDescriptor &blockInfo,
1357  Image *image)
1358  {
1359  std::list<StringContainer> imageBlock;
1360  imageBlock.push_back(filenames);
1361  CopyMetaDataToImageProperties(imageBlock, tagValueMappings_, io, blockInfo, image);
1362  }
1363 
1364  void DicomSeriesReader::CopyMetaDataToImageProperties(std::list<StringContainer> imageBlock,
1365  const gdcm::Scanner::MappingType &tagValueMappings_,
1366  DcmIoType *io,
1367  const ImageBlockDescriptor &blockInfo,
1368  Image *image)
1369  {
1370  if (!io || !image)
1371  return;
1372 
1373  StringLookupTable filesForSlices;
1374  StringLookupTable sliceLocationForSlices;
1375  StringLookupTable instanceNumberForSlices;
1376  StringLookupTable SOPInstanceNumberForSlices;
1377 
1378  auto &tagValueMappings = const_cast<gdcm::Scanner::MappingType &>(tagValueMappings_);
1379 
1380  // DICOM tags which should be added to the image properties
1381  const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location
1382 
1383  const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number
1384 
1385  const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number
1386  unsigned int timeStep(0);
1387 
1388  std::string propertyKeySliceLocation = "dicom.image.0020.1041";
1389  std::string propertyKeyInstanceNumber = "dicom.image.0020.0013";
1390  std::string propertyKeySOPInstanceNumber = "dicom.image.0008.0018";
1391 
1392  // tags for each image
1393  for (auto i = imageBlock.begin(); i != imageBlock.end(); i++, timeStep++)
1394  {
1395  const StringContainer &files = (*i);
1396  unsigned int slice(0);
1397  for (auto fIter = files.begin(); fIter != files.end(); ++fIter, ++slice)
1398  {
1399  filesForSlices.SetTableValue(slice, *fIter);
1400  gdcm::Scanner::TagToValue tagValueMapForFile = tagValueMappings[fIter->c_str()];
1401  if (tagValueMapForFile.find(tagSliceLocation) != tagValueMapForFile.end())
1402  sliceLocationForSlices.SetTableValue(slice, tagValueMapForFile[tagSliceLocation]);
1403  if (tagValueMapForFile.find(tagInstanceNumber) != tagValueMapForFile.end())
1404  instanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagInstanceNumber]);
1405  if (tagValueMapForFile.find(tagSOPInstanceNumber) != tagValueMapForFile.end())
1406  SOPInstanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagSOPInstanceNumber]);
1407  }
1408 
1409  image->SetProperty("files", StringLookupTableProperty::New(filesForSlices));
1410 
1411  // If more than one time step add postfix ".t" + timestep
1412  if (timeStep != 0)
1413  {
1414  std::ostringstream postfix;
1415  postfix << ".t" << timeStep;
1416 
1417  propertyKeySliceLocation.append(postfix.str());
1418  propertyKeyInstanceNumber.append(postfix.str());
1419  propertyKeySOPInstanceNumber.append(postfix.str());
1420  }
1421  image->SetProperty(propertyKeySliceLocation.c_str(), StringLookupTableProperty::New(sliceLocationForSlices));
1422  image->SetProperty(propertyKeyInstanceNumber.c_str(), StringLookupTableProperty::New(instanceNumberForSlices));
1423  image->SetProperty(propertyKeySOPInstanceNumber.c_str(),
1424  StringLookupTableProperty::New(SOPInstanceNumberForSlices));
1425  }
1426 
1427  // Copy tags for series, study, patient level (leave interpretation to application).
1428  // These properties will be copied to the DataNode by DicomSeriesReader.
1429 
1430  // tags for the series (we just use the one that ITK copied to its dictionary (proably that of the last slice)
1431  const itk::MetaDataDictionary &dict = io->GetMetaDataDictionary();
1433 
1434  auto dictIter = dict.Begin();
1435  while (dictIter != dict.End())
1436  {
1437  // MITK_DEBUG << "Key " << dictIter->first;
1438  std::string value;
1439  if (itk::ExposeMetaData<std::string>(dict, dictIter->first, value))
1440  {
1441  // MITK_DEBUG << "Value " << value;
1442 
1443  auto valuePosition = propertyLookup.find(dictIter->first);
1444  if (valuePosition != propertyLookup.end())
1445  {
1446  std::string propertyKey = valuePosition->second;
1447  // MITK_DEBUG << "--> " << propertyKey;
1448 
1449  image->SetProperty(propertyKey.c_str(), StringProperty::New(value));
1450  }
1451  }
1452  else
1453  {
1454  MITK_WARN << "Tag " << dictIter->first << " not read as string as expected. Ignoring...";
1455  }
1456  ++dictIter;
1457  }
1458 
1459  // copy imageblockdescriptor as properties
1460  image->SetProperty("dicomseriesreader.SOPClass", StringProperty::New(blockInfo.GetSOPClassUIDAsString()));
1461  image->SetProperty(
1462  "dicomseriesreader.ReaderImplementationLevelString",
1464  image->SetProperty("dicomseriesreader.ReaderImplementationLevel",
1466  image->SetProperty("dicomseriesreader.PixelSpacingInterpretationString",
1468  image->SetProperty("dicomseriesreader.PixelSpacingInterpretation",
1470  image->SetProperty("dicomseriesreader.MultiFrameImage", BoolProperty::New(blockInfo.IsMultiFrameImage()));
1471  image->SetProperty("dicomseriesreader.GantyTiltCorrected", BoolProperty::New(blockInfo.HasGantryTiltCorrected()));
1472  image->SetProperty("dicomseriesreader.3D+t", BoolProperty::New(blockInfo.HasMultipleTimePoints()));
1473  }
1474 
1476  {
1477  // spacing provided by ITK/GDCM
1478  Vector3D imageSpacing = image->GetGeometry()->GetSpacing();
1479  ScalarType imageSpacingX = imageSpacing[0];
1480  ScalarType imageSpacingY = imageSpacing[1];
1481 
1482  // spacing as desired by MITK (preference for "in patient", else "on detector", or "1.0/1.0")
1483  ScalarType desiredSpacingX = imageSpacingX;
1484  ScalarType desiredSpacingY = imageSpacingY;
1485  imageBlockDescriptor.GetDesiredMITKImagePixelSpacing(desiredSpacingX, desiredSpacingY);
1486 
1487  MITK_DEBUG << "Loaded spacing: " << imageSpacingX << "/" << imageSpacingY;
1488  MITK_DEBUG << "Corrected spacing: " << desiredSpacingX << "/" << desiredSpacingY;
1489 
1490  imageSpacing[0] = desiredSpacingX;
1491  imageSpacing[1] = desiredSpacingY;
1492  image->GetGeometry()->SetSpacing(imageSpacing);
1493  }
1494 
1496  DataNode &node,
1497  bool sort,
1498  bool load4D,
1499  bool correctTilt,
1500  UpdateCallBackMethod callback,
1501  Image::Pointer preLoadedImageBlock)
1502  {
1503  mitk::LocaleSwitch localeSwitch("C");
1504  std::locale previousCppLocale(std::cin.getloc());
1505  std::locale l("C");
1506  std::cin.imbue(l);
1507 
1508  ImageBlockDescriptor imageBlockDescriptor;
1509 
1510  const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image Position (Patient)
1511  const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation
1512  const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID
1513  const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP class UID
1514  const gdcm::Tag tagModality(0x0008, 0x0060); // modality
1515  const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing
1516  const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing
1517  const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames
1518 
1519  try
1520  {
1521  Image::Pointer image = preLoadedImageBlock.IsNull() ? Image::New() : preLoadedImageBlock;
1522  CallbackCommand *command = callback ? new CallbackCommand(callback) : nullptr;
1523  bool initialize_node = false;
1524 
1525  /* special case for Philips 3D+t ultrasound images */
1526  if (DicomSeriesReader::IsPhilips3DDicom(filenames.front().c_str()))
1527  {
1528  // TODO what about imageBlockDescriptor?
1529  // TODO what about preLoadedImageBlock?
1530  ReadPhilips3DDicom(filenames.front().c_str(), image);
1531  initialize_node = true;
1532  }
1533  else
1534  {
1535  /* default case: assume "normal" image blocks, possibly 3D+t */
1536  bool canLoadAs4D(true);
1537  gdcm::Scanner scanner;
1538  ScanForSliceInformation(filenames, scanner);
1539 
1540  // need non-const access for map
1541  auto &tagValueMappings = const_cast<gdcm::Scanner::MappingType &>(scanner.GetMappings());
1542 
1543  std::list<StringContainer> imageBlocks =
1544  SortIntoBlocksFor3DplusT(filenames, tagValueMappings, sort, canLoadAs4D);
1545  unsigned int volume_count = imageBlocks.size();
1546 
1547  imageBlockDescriptor.SetSeriesInstanceUID(
1548  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagSeriesInstanceUID)));
1549  imageBlockDescriptor.SetSOPClassUID(
1550  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagSOPClassUID)));
1551  imageBlockDescriptor.SetModality(
1552  DicomSeriesReader::ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagModality)));
1553  imageBlockDescriptor.SetNumberOfFrames(
1554  ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagNumberOfFrames)));
1555  imageBlockDescriptor.SetPixelSpacingInformation(
1556  ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagPixelSpacing)),
1557  ConstCharStarToString(scanner.GetValue(filenames.front().c_str(), tagImagerPixelSpacing)));
1558 
1559  GantryTiltInformation tiltInfo;
1560 
1561  // check possibility of a single slice with many timesteps. In this case, don't check for tilt, no second slice
1562  // possible
1563  if (!imageBlocks.empty() && imageBlocks.front().size() > 1 && correctTilt)
1564  {
1565  // check tiltedness here, potentially fixup ITK's loading result by shifting slice contents
1566  // check first and last position slice from tags, make some calculations to detect tilt
1567 
1568  std::string firstFilename(imageBlocks.front().front());
1569  // calculate from first and last slice to minimize rounding errors
1570  std::string secondFilename(imageBlocks.front().back());
1571 
1572  std::string imagePosition1(
1573  ConstCharStarToString(tagValueMappings[firstFilename.c_str()][tagImagePositionPatient]));
1574  std::string imageOrientation(
1575  ConstCharStarToString(tagValueMappings[firstFilename.c_str()][tagImageOrientation]));
1576  std::string imagePosition2(
1577  ConstCharStarToString(tagValueMappings[secondFilename.c_str()][tagImagePositionPatient]));
1578 
1579  bool ignoredConversionError(-42); // hard to get here, no graceful way to react
1580  Point3D origin1(DICOMStringToPoint3D(imagePosition1, ignoredConversionError));
1581  Point3D origin2(DICOMStringToPoint3D(imagePosition2, ignoredConversionError));
1582 
1583  Vector3D right;
1584  right.Fill(0.0);
1585  Vector3D up;
1586  right.Fill(0.0); // might be down as well, but it is just a name at this point
1587  DICOMStringToOrientationVectors(imageOrientation, right, up, ignoredConversionError);
1588 
1589  tiltInfo = GantryTiltInformation(origin1, origin2, right, up, filenames.size() - 1);
1590  correctTilt = tiltInfo.IsSheared() && tiltInfo.IsRegularGantryTilt();
1591  }
1592  else
1593  {
1594  correctTilt = false; // we CANNOT do that
1595  }
1596 
1597  imageBlockDescriptor.SetHasGantryTiltCorrected(correctTilt);
1598 
1599  if (volume_count == 1 || !canLoadAs4D || !load4D)
1600  {
1601  DcmIoType::Pointer io;
1602  image = MultiplexLoadDICOMByITK(
1603  imageBlocks.front(), correctTilt, tiltInfo, io, command, preLoadedImageBlock); // load first 3D block
1604 
1605  imageBlockDescriptor.AddFiles(imageBlocks.front()); // only the first part is loaded
1606  imageBlockDescriptor.SetHasMultipleTimePoints(false);
1607 
1608  FixSpacingInformation(image, imageBlockDescriptor);
1609  CopyMetaDataToImageProperties(imageBlocks.front(), scanner.GetMappings(), io, imageBlockDescriptor, image);
1610 
1611  initialize_node = true;
1612  }
1613  else if (volume_count > 1)
1614  {
1615  imageBlockDescriptor.AddFiles(filenames); // all is loaded
1616  imageBlockDescriptor.SetHasMultipleTimePoints(true);
1617 
1618  DcmIoType::Pointer io;
1619  image = MultiplexLoadDICOMByITK4D(
1620  imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command, preLoadedImageBlock);
1621 
1622  initialize_node = true;
1623  }
1624  }
1625 
1626  if (initialize_node)
1627  {
1628  // forward some image properties to node
1629  node.GetPropertyList()->ConcatenatePropertyList(image->GetPropertyList(), true);
1630 
1631  std::string patientName = "NoName";
1632  if (node.GetProperty("dicom.patient.PatientsName"))
1633  patientName = node.GetProperty("dicom.patient.PatientsName")->GetValueAsString();
1634 
1635  node.SetData(image);
1636  node.SetName(patientName);
1637  std::cin.imbue(previousCppLocale);
1638  }
1639 
1640  MITK_DEBUG << "--------------------------------------------------------------------------------";
1641  MITK_DEBUG << "DICOM files loaded (from series UID " << imageBlockDescriptor.GetSeriesInstanceUID() << "):";
1642  MITK_DEBUG << " " << imageBlockDescriptor.GetFilenames().size() << " '" << imageBlockDescriptor.GetModality()
1643  << "' files (" << imageBlockDescriptor.GetSOPClassUIDAsString() << ") loaded into 1 mitk::Image";
1644  MITK_DEBUG << " multi-frame: " << (imageBlockDescriptor.IsMultiFrameImage() ? "Yes" : "No");
1645  MITK_DEBUG << " reader support: "
1647  MITK_DEBUG << " pixel spacing type: "
1648  << PixelSpacingInterpretationToString(imageBlockDescriptor.GetPixelSpacingType()) << " "
1649  << image->GetGeometry()->GetSpacing()[0] << "/" << image->GetGeometry()->GetSpacing()[0];
1650  MITK_DEBUG << " gantry tilt corrected: " << (imageBlockDescriptor.HasGantryTiltCorrected() ? "Yes" : "No");
1651  MITK_DEBUG << " 3D+t: " << (imageBlockDescriptor.HasMultipleTimePoints() ? "Yes" : "No");
1652  MITK_DEBUG << "--------------------------------------------------------------------------------";
1653  }
1654  catch ( const std::exception &e )
1655  {
1656  // reset locale then throw up
1657  std::cin.imbue(previousCppLocale);
1658 
1659  MITK_DEBUG << "Caught exception in DicomSeriesReader::LoadDicom";
1660 
1661  throw e;
1662  }
1663  }
1664 
1665  void DicomSeriesReader::ScanForSliceInformation(const StringContainer &filenames, gdcm::Scanner &scanner)
1666  {
1667  const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image position (Patient)
1668  scanner.AddTag(tagImagePositionPatient);
1669 
1670  const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID
1671  scanner.AddTag(tagSeriesInstanceUID);
1672 
1673  const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image orientation
1674  scanner.AddTag(tagImageOrientation);
1675 
1676  const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location
1677  scanner.AddTag(tagSliceLocation);
1678 
1679  const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number
1680  scanner.AddTag(tagInstanceNumber);
1681 
1682  const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number
1683  scanner.AddTag(tagSOPInstanceNumber);
1684 
1685  const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // Pixel Spacing
1686  scanner.AddTag(tagPixelSpacing);
1687 
1688  const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // Imager Pixel Spacing
1689  scanner.AddTag(tagImagerPixelSpacing);
1690 
1691  const gdcm::Tag tagModality(0x0008, 0x0060); // Modality
1692  scanner.AddTag(tagModality);
1693 
1694  const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP Class UID
1695  scanner.AddTag(tagSOPClassUID);
1696 
1697  const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames
1698  scanner.AddTag(tagNumberOfFrames);
1699 
1700  scanner.Scan(filenames); // make available image information for each file
1701  }
1702 
1703  std::list<DicomSeriesReader::StringContainer> DicomSeriesReader::SortIntoBlocksFor3DplusT(
1704  const StringContainer &presortedFilenames,
1705  const gdcm::Scanner::MappingType &tagValueMappings,
1706  bool /*sort*/,
1707  bool &canLoadAs4D)
1708  {
1709  std::list<StringContainer> imageBlocks;
1710 
1711  // ignore sort request, because most likely re-sorting is now needed due to changes in GetSeries(bug #8022)
1712  StringContainer sorted_filenames = DicomSeriesReader::SortSeriesSlices(presortedFilenames);
1713 
1714  std::string firstPosition;
1715  unsigned int numberOfBlocks(0); // number of 3D image blocks
1716 
1717  static const gdcm::Tag tagImagePositionPatient(0x0020, 0x0032); // Image position (Patient)
1718  const gdcm::Tag tagModality(0x0008, 0x0060);
1719 
1720  // loop files to determine number of image blocks
1721  for (StringContainer::const_iterator fileIter = sorted_filenames.begin(); fileIter != sorted_filenames.end();
1722  ++fileIter)
1723  {
1724  gdcm::Scanner::TagToValue tagToValueMap = tagValueMappings.find(fileIter->c_str())->second;
1725 
1726  if (tagToValueMap.find(tagImagePositionPatient) == tagToValueMap.end())
1727  {
1728  const std::string &modality = tagToValueMap.find(tagModality)->second;
1729  if (modality.compare("RTIMAGE ") == 0 || modality.compare("RTIMAGE") == 0)
1730  {
1731  MITK_WARN << "Modality " << modality << " is not fully supported yet.";
1732  numberOfBlocks = 1;
1733  break;
1734  }
1735  else
1736  {
1737  // we expect to get images w/ missing position information ONLY as separated blocks.
1738  assert(presortedFilenames.size() == 1);
1739  numberOfBlocks = 1;
1740  break;
1741  }
1742  }
1743 
1744  std::string position = tagToValueMap.find(tagImagePositionPatient)->second;
1745  MITK_DEBUG << " " << *fileIter << " at " << position;
1746  if (firstPosition.empty())
1747  {
1748  firstPosition = position;
1749  }
1750 
1751  if (position == firstPosition)
1752  {
1753  ++numberOfBlocks;
1754  }
1755  else
1756  {
1757  break; // enough information to know the number of image blocks
1758  }
1759  }
1760 
1761  MITK_DEBUG << " ==> Assuming " << numberOfBlocks << " time steps";
1762 
1763  if (numberOfBlocks == 0)
1764  return imageBlocks; // only possible if called with no files
1765 
1766  // loop files to sort them into image blocks
1767  unsigned int numberOfExpectedSlices(0);
1768  for (unsigned int block = 0; block < numberOfBlocks; ++block)
1769  {
1770  StringContainer filesOfCurrentBlock;
1771 
1772  for (StringContainer::const_iterator fileIter = sorted_filenames.begin() + block;
1773  fileIter != sorted_filenames.end();
1774  // fileIter += numberOfBlocks) // TODO shouldn't this work? give invalid iterators on first attempts
1775  )
1776  {
1777  filesOfCurrentBlock.push_back(*fileIter);
1778  for (unsigned int b = 0; b < numberOfBlocks; ++b)
1779  {
1780  if (fileIter != sorted_filenames.end())
1781  ++fileIter;
1782  }
1783  }
1784 
1785  imageBlocks.push_back(filesOfCurrentBlock);
1786 
1787  if (block == 0)
1788  {
1789  numberOfExpectedSlices = filesOfCurrentBlock.size();
1790  }
1791  else
1792  {
1793  if (filesOfCurrentBlock.size() != numberOfExpectedSlices)
1794  {
1795  MITK_WARN
1796  << "DicomSeriesReader expected " << numberOfBlocks << " image blocks of " << numberOfExpectedSlices
1797  << " images each. Block " << block << " got " << filesOfCurrentBlock.size()
1798  << " instead. Cannot load this as 3D+t"; // TODO implement recovery (load as many slices 3D+t as much
1799  // as possible)
1800  canLoadAs4D = false;
1801  }
1802  }
1803  }
1804 
1805  return imageBlocks;
1806  }
1807 
1809  bool correctTilt,
1810  const GantryTiltInformation &tiltInfo,
1811  DcmIoType::Pointer &io,
1812  CallbackCommand *command,
1813  Image::Pointer preLoadedImageBlock)
1814  {
1815  io = DcmIoType::New();
1816  io->SetFileName(filenames.front().c_str());
1817  io->ReadImageInformation();
1818 
1819  if (io->GetPixelType() == itk::ImageIOBase::SCALAR)
1820  {
1821  return MultiplexLoadDICOMByITKScalar(filenames, correctTilt, tiltInfo, io, command, preLoadedImageBlock);
1822  }
1823  else if (io->GetPixelType() == itk::ImageIOBase::RGB)
1824  {
1825  return MultiplexLoadDICOMByITKRGBPixel(filenames, correctTilt, tiltInfo, io, command, preLoadedImageBlock);
1826  }
1827  else
1828  {
1829  return nullptr;
1830  }
1831  }
1832 
1833  Image::Pointer DicomSeriesReader::MultiplexLoadDICOMByITK4D(std::list<StringContainer> &imageBlocks,
1834  ImageBlockDescriptor imageBlockDescriptor,
1835  bool correctTilt,
1836  const GantryTiltInformation &tiltInfo,
1837  DcmIoType::Pointer &io,
1838  CallbackCommand *command,
1839  Image::Pointer preLoadedImageBlock)
1840  {
1841  io = DcmIoType::New();
1842  io->SetFileName(imageBlocks.front().front().c_str());
1843  io->ReadImageInformation();
1844 
1845  if (io->GetPixelType() == itk::ImageIOBase::SCALAR)
1846  {
1848  imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command, preLoadedImageBlock);
1849  }
1850  else if (io->GetPixelType() == itk::ImageIOBase::RGB)
1851  {
1853  imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command, preLoadedImageBlock);
1854  }
1855  else
1856  {
1857  return nullptr;
1858  }
1859  }
1860 } // end namespace mitk
PixelSpacingInterpretation
How the mitk::Image spacing should be interpreted.
mitk::BaseProperty * GetProperty(const char *propertyKey, const mitk::BaseRenderer *renderer=nullptr, bool fallBackOnDataProperties=true) const
Get the property (instance of BaseProperty) with key propertyKey from the PropertyList of the rendere...
StringContainer GetBlockFilenames()
Grouping result, all same origin-to-origin distance w/o gaps.
std::map< std::string, std::string > TagToPropertyMapType
Maps DICOM tags to MITK properties.
static itk::SmartPointer< Image > MultiplexLoadDICOMByITK4DScalar(std::list< StringContainer > &imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
static bool DICOMStringToSpacing(const std::string &s, ScalarType &spacingX, ScalarType &spacingY)
Safely convert a string into pixel spacing x and y.
bool IsSheared() const
Whether the slices were sheared.
void SetSpacing(const mitk::Vector3D &aSpacing, bool enforceSetSpacing=false)
Set the spacing (m_Spacing).
PixelSpacingInterpretation GetPixelSpacingType() const
How the mitk::Image spacing can meaningfully be interpreted.
Return type of GetSeries, describes a logical group of files.
void SetTableValue(IdentifierType id, ValueType value)
static Point3D DICOMStringToPoint3D(const std::string &s, bool &successful)
Convert DICOM string describing a point to Point3D.
static void ScanForSliceInformation(const StringContainer &filenames, gdcm::Scanner &scanner)
Scan for slice image information.
itk::Image< unsigned char, 3 > ImageType
#define MITK_ERROR
Definition: mitkLogMacros.h:20
double ScalarType
static void FixSpacingInformation(Image *image, const ImageBlockDescriptor &imageBlockDescriptor)
void UndoPrematureGrouping()
Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption.
static std::string GetConfigurationString()
Provide combination of preprocessor defines that was active during compilation.
virtual void SetData(mitk::BaseData *baseData)
Set the data object (instance of BaseData, e.g., an Image) managed by this DataNode.
std::vector< std::string > StringContainer
Lists of filenames.
static bool ReadPhilips3DDicom(const std::string &filename, itk::SmartPointer< Image > output_image)
Read a Philips3D ultrasound DICOM file and put into an mitk::Image.
#define MITK_DEBUG
Definition: mitkLogMacros.h:22
static DataNode::Pointer LoadDicomSeries(const StringContainer &filenames, bool sort=true, bool load4D=true, bool correctGantryTilt=true, UpdateCallBackMethod callback=nullptr, itk::SmartPointer< Image > preLoadedImageBlock=nullptr)
DataCollection - Class to facilitate loading/accessing structured data.
void ConcatenatePropertyList(PropertyList *pList, bool replace=false)
Set a property object in the list/map by reference.
std::string GetSOPClassUIDAsString() const
SOP Class UID as readable string (Computed Tomography Image Storage, Secondary Capture Image Storage...
static itk::SmartPointer< Image > MultiplexLoadDICOMByITKRGBPixel(const StringContainer &, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
void(* UpdateCallBackMethod)(float)
Interface for the progress callback.
StringContainer GetFilenames() const
List of files in this group.
StringContainer GetUnsortedFilenames()
Remaining files, which could not be grouped.
static std::string CreateSeriesIdentifierPart(gdcm::Scanner::TagToValue &tagValueMap, const gdcm::Tag &tag)
Helper for CreateMoreUniqueSeriesIdentifier.
static std::string ReaderImplementationLevelToString(const ReaderImplementationLevel &enumValue)
static itk::SmartPointer< Image > MultiplexLoadDICOMByITK4DRGBPixel(std::list< StringContainer > &imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
static std::string ConstCharStarToString(const char *s)
Safely convert const char* to std::string.
loader code is implemented but not accompanied by tests
void SetProperty(const std::string &propertyKey, BaseProperty *property, const std::string &contextName="", bool fallBackOnDefaultContext=false) override
Add new or change existent property.
std::string GetSeriesInstanceUID() const
The Series Instance UID.
static Pointer New()
static std::string PixelSpacingInterpretationToString(const PixelSpacingInterpretation &enumValue)
static std::string IDifyTagValue(const std::string &value)
Helper for CreateMoreUniqueSeriesIdentifier.
BaseData * GetData() const
Get the data object (instance of BaseData, e.g., an Image) managed by this DataNode.
virtual std::string GetValueAsString() const
#define MITK_WARN
Definition: mitkLogMacros.h:19
static void LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback, itk::SmartPointer< Image > preLoadedImageBlock)
Performs actual loading of a series and creates an image having the specified pixel type...
Convenience class to temporarily change the current locale.
static std::list< StringContainer > SortIntoBlocksFor3DplusT(const StringContainer &presortedFilenames, const gdcm::Scanner::MappingType &tagValueMappings_, bool sort, bool &canLoadAs4D)
Sort files into time step blocks of a 3D+t image.
static Pointer New()
ReaderImplementationLevel GetReaderImplementationLevel() const
Confidence of the reader that this block can be read successfully.
#define mitkThrow()
static FileNamesGrouping GetSeries(const std::string &dir, bool groupImagesWithGantryTilt, const StringContainer &restrictions=StringContainer())
see other GetSeries().
static void DICOMStringToOrientationVectors(const std::string &s, Vector3D &right, Vector3D &up, bool &successful)
Convert DICOM string describing a point two Vector3D.
void AddFileToSortedBlock(const std::string &filename)
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.
Image class for storing images.
Definition: mitkImage.h:72
virtual void SetName(const char *name)
Extra convenience access method to set the name of an object.
Definition: mitkDataNode.h:403
bool ContainsGantryTilt()
Wheter or not the grouped result contain a gantry tilt.
Progress callback for DicomSeriesReader.
std::string GetModality() const
Series Modality (CT, MR, etc.)
mitk::Image::Pointer image
static bool IsPhilips3DDicom(const std::string &filename)
Checks if a specific file is a Philips3D ultrasound DICOM file.
static itk::SmartPointer< Image > MultiplexLoadDICOMByITKScalar(const StringContainer &, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
static bool IsDicom(const std::string &filename)
Checks if a specific file contains DICOM data.
static Pointer New()
std::string GetImageBlockUID() const
A unique ID describing this bloc (enhanced Series Instance UID).
mitk::PropertyList * GetPropertyList(const mitk::BaseRenderer *renderer=nullptr) const
Get the PropertyList of the renderer. If renderer is nullptr, the BaseRenderer-independent PropertyLi...
Return type of DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption.
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
static void CopyMetaDataToImageProperties(StringContainer filenames, const gdcm::Scanner::MappingType &tagValueMappings_, DcmIoType *io, const ImageBlockDescriptor &blockInfo, Image *image)
Copy information about files and DICOM tags from ITK&#39;s MetaDataDictionary and from the list of input ...
static itk::SmartPointer< Image > MultiplexLoadDICOMByITK4D(std::list< StringContainer > &imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
void FlagGantryTilt()
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.
MITKCORE_EXPORT const ScalarType eps
static const TagToPropertyMapType & GetDICOMTagsToMITKPropertyMap()
Map between DICOM tags and MITK properties.
bool IsRegularGantryTilt() const
Whether the shearing is a gantry tilt or more complicated.
void AddFileToUnsortedBlock(const std::string &filename)
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.
std::map< std::string, ImageBlockDescriptor > FileNamesGrouping
bool IsMultiFrameImage() const
Multi-frame image(s) or not.
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:143
static std::string CreateMoreUniqueSeriesIdentifier(gdcm::Scanner::TagToValue &tagValueMap)
Construct a UID that takes into account sorting criteria from GetSeries().
static Pointer New()
static itk::SmartPointer< Image > MultiplexLoadDICOMByITK(const StringContainer &, bool correctTilt, const GantryTiltInformation &tiltInfo, DcmIoType::Pointer &io, CallbackCommand *command, itk::SmartPointer< Image > preLoadedImageBlock)
static StringContainer SortSeriesSlices(const StringContainer &unsortedFilenames)
Sort a set of file names in an order that is meaningful for loading them into an mitk::Image.
ReaderImplementationLevel
Describes how well the reader is tested for a certain file type.
static bool GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2)
Defines spatial sorting for sorting by GDCM 2.
Class for nodes of the DataTree.
Definition: mitkDataNode.h:57
bool HasGantryTiltCorrected() const
Whether or not the block contains a gantry tilt which will be "corrected" during loading.
static SliceGroupingAnalysisResult AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer &files, bool groupsOfSimilarImages, const gdcm::Scanner::MappingType &tagValueMappings_)
Ensure an equal z-spacing for a group of files.