Medical Imaging Interaction Toolkit  2018.4.99-b7f3afaa
Medical Imaging Interaction Toolkit
mitkEquiDistantBlocksSorter.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 //#define MBILOG_ENABLE_DEBUG
14 
16 
19 {
20 }
21 
25 {
26  return m_GroupedFiles;
27 }
28 
32 {
33  return m_UnsortedFiles;
34 }
35 
36 bool
39 {
40  return m_TiltInfo.IsRegularGantryTilt();
41 }
42 
43 void
46 {
47  m_GroupedFiles.push_back( dataset );
48 }
49 
50 void
53 {
54  m_UnsortedFiles.push_back( dataset );
55 }
56 
57 void
60 {
61  m_UnsortedFiles.insert( m_UnsortedFiles.end(), datasets.begin(), datasets.end() );
62 }
63 
64 void
66 ::SetFirstFilenameOfBlock(const std::string& filename)
67 {
68  m_FirstFilenameOfBlock = filename;
69 }
70 
71 std::string
74 {
75  return m_FirstFilenameOfBlock;
76 }
77 
78 void
80 ::SetLastFilenameOfBlock(const std::string& filename)
81 {
82  m_LastFilenameOfBlock = filename;
83 }
84 
85 std::string
88 {
89  return m_LastFilenameOfBlock;
90 }
91 
92 
93 void
96 {
97  m_TiltInfo = tiltInfo;
98 }
99 
103 {
104  return m_TiltInfo;
105 }
106 
107 void
110 {
111  assert( !m_GroupedFiles.empty() );
112  m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() );
113  m_GroupedFiles.pop_back();
114  m_TiltInfo = GantryTiltInformation();
115 }
116 
117 // ------------------------ end helper class
118 
122 ,m_AcceptTilt(false)
126 {
127 }
128 
131 :DICOMDatasetSorter(other)
136 {
137 }
138 
141 {
142 }
143 
144 bool
146 ::operator==(const DICOMDatasetSorter& other) const
147 {
148  if (const auto* otherSelf = dynamic_cast<const EquiDistantBlocksSorter*>(&other))
149  {
150  return this->m_AcceptTilt == otherSelf->m_AcceptTilt
151  && this->m_ToleratedOriginOffsetIsAbsolute == otherSelf->m_ToleratedOriginOffsetIsAbsolute
152  && this->m_AcceptTwoSlicesGroups == otherSelf->m_AcceptTwoSlicesGroups
153  && (fabs(this->m_ToleratedOriginOffset - otherSelf->m_ToleratedOriginOffset) < eps);
154  }
155  else
156  {
157  return false;
158  }
159 }
160 
161 void
163 ::PrintConfiguration(std::ostream& os, const std::string& indent) const
164 {
165  std::stringstream ts;
167  {
168  ts << "adaptive";
169  }
170  else
171  {
172  ts << m_ToleratedOriginOffset << "mm";
173  }
174 
175  os << indent << "Sort into blocks of equidistant, well-aligned (tolerance "
176  << ts.str() << ") slices "
177  << (m_AcceptTilt ? "(accepting a gantry tilt)" : "")
178  << std::endl;
179 }
180 
181 
182 void
184 ::SetAcceptTilt(bool accept)
185 {
186  m_AcceptTilt = accept;
187 }
188 
189 
190 bool
193 {
194  return m_AcceptTilt;
195 }
196 
197 void
200 {
201  m_AcceptTwoSlicesGroups = accept;
202 }
203 
204 bool
207 {
209 }
210 
211 
215 {
216  if (this != &other)
217  {
219  m_AcceptTilt = other.m_AcceptTilt;
223  }
224  return *this;
225 }
226 
230 {
231  DICOMTagList tags;
232  tags.push_back( DICOMTag(0x0020, 0x0032) ); // ImagePositionPatient
233  tags.push_back( DICOMTag(0x0020, 0x0037) ); // ImageOrientationPatient
234  tags.push_back( DICOMTag(0x0018, 0x1120) ); // GantryDetectorTilt
235 
236  return tags;
237 }
238 
239 void
242 {
243  DICOMDatasetList remainingInput = GetInput(); // copy
244 
245  typedef std::list<DICOMDatasetList> OutputListType;
246  OutputListType outputs;
247 
248  m_SliceGroupingResults.clear();
249 
250  while (!remainingInput.empty()) // repeat until all files are grouped somehow
251  {
253 
254 #ifdef MBILOG_ENABLE_DEBUG
255  DICOMDatasetList inBlock = regularBlock.GetBlockDatasets();
256  DICOMDatasetList laterBlock = regularBlock.GetUnsortedDatasets();
257  MITK_DEBUG << "Result: sorted 3D group with " << inBlock.size() << " files";
258  for (DICOMDatasetList::const_iterator diter = inBlock.cbegin(); diter != inBlock.cend(); ++diter)
259  MITK_DEBUG << " IN " << (*diter)->GetFilenameIfAvailable();
260  for (DICOMDatasetList::const_iterator diter = laterBlock.cbegin(); diter != laterBlock.cend(); ++diter)
261  MITK_DEBUG << " OUT " << (*diter)->GetFilenameIfAvailable();
262 #endif // MBILOG_ENABLE_DEBUG
263 
264  outputs.push_back( regularBlock.GetBlockDatasets() );
265  m_SliceGroupingResults.push_back( regularBlock );
266  remainingInput = regularBlock.GetUnsortedDatasets();
267  }
268 
269  unsigned int numberOfOutputs = outputs.size();
270  this->SetNumberOfOutputs(numberOfOutputs);
271 
272  unsigned int outputIndex(0);
273  for (auto oIter = outputs.cbegin();
274  oIter != outputs.cend();
275  ++outputIndex, ++oIter)
276  {
277  this->SetOutput(outputIndex, *oIter);
278  }
279 }
280 
281 void
283 ::SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistance)
284 {
285  m_ToleratedOriginOffset = fractionOfInterSliceDistance;
287 
288  if (m_ToleratedOriginOffset < 0.0)
289  {
290  MITK_WARN << "Call SetToleratedOriginOffsetToAdaptive() only with positive numbers between 0.0 and 1.0, read documentation!";
291  }
292 
293  if (m_ToleratedOriginOffset > 0.5)
294  {
295  MITK_WARN << "EquiDistantBlocksSorter is now accepting large errors, take care of measurements, they could appear at unprecise locations!";
296  }
297 }
298 
299 void
301 ::SetToleratedOriginOffset(double millimeters)
302 {
303  m_ToleratedOriginOffset = millimeters;
305  if (m_ToleratedOriginOffset < 0.0)
306  {
307  MITK_WARN << "Negative tolerance set to SetToleratedOriginOffset()!";
308  }
309 }
310 
311 double
314 {
316 }
317 
318 bool
321 {
323 }
324 
325 
326 std::string
329 {
330  return s ? std::string(s) : std::string();
331 }
332 
336  const DICOMDatasetList& datasets,
337  bool groupImagesWithGantryTilt)
338 {
339  // result.first = files that fit ITK's assumption
340  // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again
342 
343  const DICOMTag tagImagePositionPatient = DICOMTag(0x0020,0x0032); // Image Position (Patient)
344  const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation
345 
346  Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0);
347  bool fromFirstToSecondOriginInitialized(false);
348  Point3D thisOrigin;
349  thisOrigin.Fill(0.0f);
350  Point3D lastOrigin;
351  lastOrigin.Fill(0.0f);
352  Point3D lastDifferentOrigin;
353  lastDifferentOrigin.Fill(0.0f);
354 
355  bool lastOriginInitialized(false);
356 
357  MITK_DEBUG << "--------------------------------------------------------------------------------";
358  MITK_DEBUG << "Analyzing " << datasets.size() << " files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")";
359  unsigned int fileIndex(0);
360  double toleratedOriginError(0.005); // default: max. 1/10mm error when measurement crosses 20 slices in z direction (too strict? we don't know better)
361  for (auto dsIter = datasets.cbegin();
362  dsIter != datasets.cend();
363  ++dsIter, ++fileIndex)
364  {
365  bool fileFitsIntoPattern(false);
366  std::string thisOriginString;
367  // Read tag value into point3D. PLEASE replace this by appropriate GDCM code if you figure out how to do that
368  thisOriginString = (*dsIter)->GetTagValueAsString(tagImagePositionPatient).value;
369 
370  if (thisOriginString.empty())
371  {
372  // don't let such files be in a common group. Everything without position information will be loaded as a single slice:
373  // with standard DICOM files this can happen to: CR, DX, SC
374  MITK_DEBUG << " ==> Sort away " << *dsIter << " for later analysis (no position information)"; // we already have one occupying this position
375 
376  if ( result.GetBlockDatasets().empty() ) // nothing WITH position information yet
377  {
378  // ==> this is a group of its own, stop processing, come back later
379  result.AddFileToSortedBlock( *dsIter );
380 
381  DICOMDatasetList remainingFiles;
382  remainingFiles.insert( remainingFiles.end(), dsIter+1, datasets.end() );
383  result.AddFilesToUnsortedBlock( remainingFiles );
384 
385  fileFitsIntoPattern = false;
386  break; // no files anymore
387  }
388  else
389  {
390  // ==> this does not match, consider later
391  result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis
392  fileFitsIntoPattern = false;
393  continue; // next file
394  }
395  }
396 
397  bool ignoredConversionError(-42); // hard to get here, no graceful way to react
398  thisOrigin = DICOMStringToPoint3D( thisOriginString, ignoredConversionError );
399 
400  MITK_DEBUG << " " << fileIndex << " " << (*dsIter)->GetFilenameIfAvailable()
401  << " at "
402  /* << thisOriginString */ << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")";
403 
404  if ( lastOriginInitialized && (thisOrigin == lastOrigin) )
405  {
406  MITK_DEBUG << " ==> Sort away " << *dsIter << " for separate time step"; // we already have one occupying this position
407  result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis
408  fileFitsIntoPattern = false;
409  }
410  else
411  {
412  if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position
413  {
414  fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin;
415  fromFirstToSecondOriginInitialized = true;
416 
417  // classic mode without tolerance!
419  {
420  MITK_DEBUG << "Distance of two slices: " << fromFirstToSecondOrigin.GetNorm() << "mm";
421  toleratedOriginError =
422  fromFirstToSecondOrigin.GetNorm() * 0.3; // a third of the slice distance
423  // (less than half, which would mean that a slice is displayed where another slice should actually be)
424  }
425  else
426  {
427  toleratedOriginError = m_ToleratedOriginOffset;
428  }
429  MITK_DEBUG << "Accepting errors in actual versus expected origin up to " << toleratedOriginError << "mm";
430 
431  // Here we calculate if this slice and the previous one are well aligned,
432  // i.e. we test if the previous origin is on a line through the current
433  // origin, directed into the normal direction of the current slice.
434 
435  // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry,
436  // which cannot be simply loaded into a single mitk::Image at the moment.
437  // For this case, we flag this finding in the result and DicomSeriesReader
438  // can correct for that later.
439 
440  Vector3D right; right.Fill(0.0);
441  Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point
442  std::string orientationValue = (*dsIter)->GetTagValueAsString( tagImageOrientation ).value;
443  DICOMStringToOrientationVectors( orientationValue, right, up, ignoredConversionError );
444 
445  GantryTiltInformation tiltInfo( lastDifferentOrigin, thisOrigin, right, up, 1 );
446 
447  if ( tiltInfo.IsSheared() )
448  {
449  /* optimistic approach, accepting gantry tilt: save file for later, check all further files */
450 
451  // at this point we have TWO slices analyzed! if they are the only two files, we still split, because there is no third to verify our tilting assumption.
452  // later with a third being available, we must check if the initial tilting vector is still valid. if yes, continue.
453  // if NO, we need to split the already sorted part (result.first) and the currently analyzed file (*dsIter)
454 
455  // tell apart gantry tilt from overall skewedness
456  // sort out irregularly sheared slices, that IS NOT tilting
457 
458  if ( groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt() )
459  {
460  assert(!datasets.empty());
461 
462  result.FlagGantryTilt(tiltInfo);
463  result.AddFileToSortedBlock( *dsIter ); // this file is good for current block
464  result.SetFirstFilenameOfBlock( datasets.front()->GetFilenameIfAvailable() );
465  result.SetLastFilenameOfBlock( datasets.back()->GetFilenameIfAvailable() );
466  fileFitsIntoPattern = true;
467  }
468  else // caller does not want tilt compensation OR shearing is more complicated than tilt
469  {
470  result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis
471  fileFitsIntoPattern = false;
472  }
473  }
474  else // not sheared
475  {
476  result.AddFileToSortedBlock( *dsIter ); // this file is good for current block
477  fileFitsIntoPattern = true;
478  }
479  }
480  else if (fromFirstToSecondOriginInitialized) // we already know the offset between slices
481  {
482  Point3D assumedOrigin = lastDifferentOrigin + fromFirstToSecondOrigin;
483 
484  Vector3D originError = assumedOrigin - thisOrigin;
485  double norm = originError.GetNorm();
486 
487  if (norm > toleratedOriginError)
488  {
489  MITK_DEBUG << " File does not fit into the inter-slice distance pattern (diff = "
490  << norm << ", allowed "
491  << toleratedOriginError << ").";
492  MITK_DEBUG << " Expected position (" << assumedOrigin[0] << ","
493  << assumedOrigin[1] << ","
494  << assumedOrigin[2] << "), got position ("
495  << thisOrigin[0] << ","
496  << thisOrigin[1] << ","
497  << thisOrigin[2] << ")";
498  MITK_DEBUG << " ==> Sort away " << *dsIter << " for later analysis";
499 
500  // At this point we know we deviated from the expectation of ITK's ImageSeriesReader
501  // We split the input file list at this point, i.e. all files up to this one (excluding it)
502  // are returned as group 1, the remaining files (including the faulty one) are group 2
503 
504  /* Optimistic approach: check if any of the remaining slices fits in */
505  result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis
506  fileFitsIntoPattern = false;
507  }
508  else
509  {
510  result.AddFileToSortedBlock( *dsIter ); // this file is good for current block
511  fileFitsIntoPattern = true;
512  }
513  }
514  else // this should be the very first slice
515  {
516  result.AddFileToSortedBlock( *dsIter ); // this file is good for current block
517  fileFitsIntoPattern = true;
518  }
519  }
520 
521  // record current origin for reference in later iterations
522  if ( !lastOriginInitialized || ( fileFitsIntoPattern && (thisOrigin != lastOrigin) ) )
523  {
524  lastDifferentOrigin = thisOrigin;
525  }
526 
527  lastOrigin = thisOrigin;
528  lastOriginInitialized = true;
529  }
530 
531  if ( result.ContainsGantryTilt() )
532  {
533  // check here how many files were grouped.
534  // IF it was only two files AND we assume tiltedness (e.g. save "distance")
535  // THEN we would want to also split the two previous files (simple) because
536  // we don't have any reason to assume they belong together
537 
538  // Above behavior can be configured via m_AcceptTwoSlicesGroups, the default being "do accept"
539  if ( result.GetBlockDatasets().size() == 2 && !m_AcceptTwoSlicesGroups )
540  {
541  result.UndoPrematureGrouping();
542  }
543  }
544 
545  // update tilt info to get maximum precision
546  // earlier, tilt was only calculated from first and second slice.
547  // now that we know the whole range, we can re-calculate using the very first and last slice
548  if ( result.ContainsGantryTilt() && result.GetBlockDatasets().size() > 1 )
549  {
550  try
551  {
552  DICOMDatasetList datasets = result.GetBlockDatasets();
553  DICOMDatasetAccess* firstDataset = datasets.front();
554  DICOMDatasetAccess* lastDataset = datasets.back();
555  unsigned int numberOfSlicesApart = datasets.size() - 1;
556 
557  std::string orientationString = firstDataset->GetTagValueAsString( tagImageOrientation ).value;
558  std::string firstOriginString = firstDataset->GetTagValueAsString( tagImagePositionPatient ).value;
559  std::string lastOriginString = lastDataset->GetTagValueAsString( tagImagePositionPatient ).value;
560 
561  result.FlagGantryTilt( GantryTiltInformation::MakeFromTagValues( firstOriginString, lastOriginString, orientationString, numberOfSlicesApart ));
562  }
563  catch (...)
564  {
565  // just do not flag anything, we are ok
566  }
567  }
568 
569  return result;
570 }
SliceGroupingAnalysisResult AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const DICOMDatasetList &files, bool groupsOfSimilarImages)
Ensure an equal z-spacing for a group of files.
bool ContainsGantryTilt()
Wheter or not the grouped result contain a gantry tilt.
EquiDistantBlocksSorter & operator=(const EquiDistantBlocksSorter &other)
void SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistanct=0.3)
See class description and SetToleratedOriginOffset().
DICOMDatasetList GetUnsortedDatasets()
Remaining files, which could not be grouped.
static GantryTiltInformation MakeFromTagValues(const std::string &origin1String, const std::string &origin2String, const std::string &orientationString, unsigned int numberOfSlicesApart)
Factory method to create GantryTiltInformation from tag values (strings).
The sorting/splitting building-block of DICOMITKSeriesGDCMReader.
void Sort() override
Delegates work to AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). AnalyzeFileForITKImageSeries...
void SetAcceptTilt(bool accept)
Whether or not to accept images from a tilted acquisition in a single output group.
Split inputs into blocks of equidistant slices (for use in DICOMITKSeriesGDCMReader).
std::vector< DICOMTag > DICOMTagList
Definition: mitkDICOMTag.h:59
void DICOMStringToOrientationVectors(const std::string &s, Vector3D &right, Vector3D &up, bool &successful)
Convert DICOM string describing a point two Vector3D.
Representation of a DICOM tag.
Definition: mitkDICOMTag.h:32
#define MITK_DEBUG
Definition: mitkLogMacros.h:22
DICOMTagList GetTagsOfInterest() override
Return the tags of interest (to facilitate scanning)
void SetNumberOfOutputs(unsigned int numberOfOutputs)
std::string ConstCharStarToString(const char *s)
Safely convert const char* to std::string.
const DICOMDatasetList & GetInput() const
Input for sorting.
Interface to datasets that is presented to sorting classes such as DICOMDatasetSorter.
Return type of AnalyzeFileForITKImageSeriesReaderSpacingAssumption().
const GantryTiltInformation & GetTiltInfo() const
Detailed description of gantry tilt.
DICOMDatasetList GetBlockDatasets()
Grouping result, all same origin-to-origin distance w/o gaps.
#define MITK_WARN
Definition: mitkLogMacros.h:19
virtual DICOMDatasetFinding GetTagValueAsString(const DICOMTag &tag) const =0
Return a DICOMDatasetFinding instance of the tag. The return containes (if valid) the raw value of th...
void AddFileToSortedBlock(DICOMDatasetAccess *dataset)
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.
Gantry tilt analysis result.
bool IsRegularGantryTilt() const
Whether the shearing is a gantry tilt or more complicated.
MITKCORE_EXPORT const ScalarType eps
Point3D DICOMStringToPoint3D(const std::string &s, bool &successful)
Convert DICOM string describing a point to Point3D.
void SetOutput(unsigned int index, const DICOMDatasetList &output)
DICOMDatasetSorter & operator=(const DICOMDatasetSorter &other)
void UndoPrematureGrouping()
Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption.
std::vector< DICOMDatasetAccess * > DICOMDatasetList
bool operator==(const DICOMDatasetSorter &other) const override
bool IsSheared() const
Whether the slices were sheared.
void SetToleratedOriginOffset(double millimeters=0.005)
See class description and SetToleratedOriginOffsetToAdaptive().
void PrintConfiguration(std::ostream &os, const std::string &indent="") const override
Print configuration details into stream.
void FlagGantryTilt(const GantryTiltInformation &tiltInfo)
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.
void AddFileToUnsortedBlock(DICOMDatasetAccess *dataset)
Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only.