Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkSegmentationInterpolationController.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,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
18 
19 #include "mitkImageCast.h"
20 #include "mitkImageReadAccessor.h"
21 #include "mitkImageTimeSelector.h"
22 #include <mitkExtractSliceFilter.h>
23 #include <mitkImageAccessByItk.h>
24 //#include <mitkPlaneGeometry.h>
25 
27 
28 #include <itkCommand.h>
29 #include <itkImage.h>
30 #include <itkImageSliceConstIteratorWithIndex.h>
31 
34 
36  const Image *image)
37 {
38  auto iter = s_InterpolatorForImage.find(image);
39  if (iter != s_InterpolatorForImage.end())
40  {
41  return iter->second;
42  }
43  else
44  {
45  return nullptr;
46  }
47 }
48 
50 {
51 }
52 
54 {
55  m_2DInterpolationActivated = status;
56 }
57 
59 {
61 
62  if (m_Instance.IsNull())
63  {
65  }
66  return m_Instance;
67 }
68 
70 {
71  // remove this from the list of interpolators
72  for (auto iter = s_InterpolatorForImage.begin(); iter != s_InterpolatorForImage.end(); ++iter)
73  {
74  if (iter->second == this)
75  {
76  s_InterpolatorForImage.erase(iter);
77  break;
78  }
79  }
80 }
81 
83 {
84  if (!m_BlockModified && m_Segmentation.IsNotNull() && m_2DInterpolationActivated)
85  {
86  SetSegmentationVolume(m_Segmentation);
87  }
88 }
89 
91 {
92  m_BlockModified = block;
93 }
94 
96 {
97  // clear old information (remove all time steps
98  m_SegmentationCountInSlice.clear();
99 
100  // delete this from the list of interpolators
101  auto iter = s_InterpolatorForImage.find(segmentation);
102  if (iter != s_InterpolatorForImage.end())
103  {
104  s_InterpolatorForImage.erase(iter);
105  }
106 
107  if (!segmentation)
108  return;
109  if (segmentation->GetDimension() > 4 || segmentation->GetDimension() < 3)
110  {
111  itkExceptionMacro("SegmentationInterpolationController needs a 3D-segmentation or 3D+t, not 2D.");
112  }
113 
114  if (m_Segmentation != segmentation)
115  {
116  // observe Modified() event of image
119  command->SetCallbackFunction(this, &SegmentationInterpolationController::OnImageModified);
120  segmentation->AddObserver(itk::ModifiedEvent(), command);
121  }
122 
123  m_Segmentation = segmentation;
124 
125  m_SegmentationCountInSlice.resize(m_Segmentation->GetTimeSteps());
126  for (unsigned int timeStep = 0; timeStep < m_Segmentation->GetTimeSteps(); ++timeStep)
127  {
128  m_SegmentationCountInSlice[timeStep].resize(3);
129  for (unsigned int dim = 0; dim < 3; ++dim)
130  {
131  m_SegmentationCountInSlice[timeStep][dim].clear();
132  m_SegmentationCountInSlice[timeStep][dim].resize(m_Segmentation->GetDimension(dim));
133  m_SegmentationCountInSlice[timeStep][dim].assign(m_Segmentation->GetDimension(dim), 0);
134  }
135  }
136 
137  s_InterpolatorForImage.insert(std::make_pair(m_Segmentation, this));
138 
139  // for all timesteps
140  // scan whole image
141  for (unsigned int timeStep = 0; timeStep < m_Segmentation->GetTimeSteps(); ++timeStep)
142  {
144  timeSelector->SetInput(m_Segmentation);
145  timeSelector->SetTimeNr(timeStep);
146  timeSelector->UpdateLargestPossibleRegion();
147  Image::Pointer segmentation3D = timeSelector->GetOutput();
148  AccessFixedDimensionByItk_2(segmentation3D, ScanWholeVolume, 3, m_Segmentation, timeStep);
149  }
150 
151  // PrintStatus();
152 
153  SetReferenceVolume(m_ReferenceImage);
154 
155  Modified();
156 }
157 
159 {
160  m_ReferenceImage = referenceImage;
161 
162  if (m_ReferenceImage.IsNull())
163  return; // no image set - ignore it then
164  assert(m_Segmentation.IsNotNull()); // should never happen
165 
166  // ensure the reference image has the same dimensionality and extents as the segmentation image
167  if (m_ReferenceImage.IsNull() || m_Segmentation.IsNull() ||
168  m_ReferenceImage->GetDimension() != m_Segmentation->GetDimension() ||
169  m_ReferenceImage->GetPixelType().GetNumberOfComponents() != 1 ||
170  m_Segmentation->GetPixelType().GetNumberOfComponents() != 1)
171  {
172  MITK_WARN << "Segmentation image has different image characteristics than reference image." << std::endl;
173  m_ReferenceImage = nullptr;
174  return;
175  }
176 
177  for (unsigned int dim = 0; dim < m_Segmentation->GetDimension(); ++dim)
178  if (m_ReferenceImage->GetDimension(dim) != m_Segmentation->GetDimension(dim))
179  {
180  MITK_WARN << "original patient image does not match segmentation (different extent in dimension " << dim
181  << "), ignoring patient image" << std::endl;
182  m_ReferenceImage = nullptr;
183  return;
184  }
185 }
186 
187 void mitk::SegmentationInterpolationController::SetChangedVolume(const Image *sliceDiff, unsigned int timeStep)
188 {
189  if (!sliceDiff)
190  return;
191  if (sliceDiff->GetDimension() != 3)
192  return;
193 
194  AccessFixedDimensionByItk_1(sliceDiff, ScanChangedVolume, 3, timeStep);
195 
196  // PrintStatus();
197  Modified();
198 }
199 
201  unsigned int sliceDimension,
202  unsigned int sliceIndex,
203  unsigned int timeStep)
204 {
205  if (!sliceDiff)
206  return;
207  if (sliceDimension > 2)
208  return;
209  if (timeStep >= m_SegmentationCountInSlice.size())
210  return;
211  if (sliceIndex >= m_SegmentationCountInSlice[timeStep][sliceDimension].size())
212  return;
213 
214  unsigned int dim0(0);
215  unsigned int dim1(1);
216 
217  // determine the other two dimensions
218  switch (sliceDimension)
219  {
220  default:
221  case 2:
222  dim0 = 0;
223  dim1 = 1;
224  break;
225  case 1:
226  dim0 = 0;
227  dim1 = 2;
228  break;
229  case 0:
230  dim0 = 1;
231  dim1 = 2;
232  break;
233  }
234 
235  // mitkIpPicDescriptor* rawSlice = const_cast<Image*>(sliceDiff)->GetSliceData()->GetPicDescriptor(); // we promise
236  // not to change anything!
237 
238  mitk::ImageReadAccessor readAccess(sliceDiff);
239  unsigned char *rawSlice = (unsigned char *)readAccess.GetData();
240  if (!rawSlice)
241  return;
242 
244  sliceDiff, ScanChangedSlice, 2, SetChangedSliceOptions(sliceDimension, sliceIndex, dim0, dim1, timeStep, rawSlice));
245 
246  // PrintStatus();
247 
248  Modified();
249 }
250 
251 template <typename DATATYPE>
253  const SetChangedSliceOptions &options)
254 {
255  DATATYPE *pixelData((DATATYPE *)options.pixelData);
256 
257  unsigned int timeStep(options.timeStep);
258 
259  unsigned int sliceDimension(options.sliceDimension);
260  unsigned int sliceIndex(options.sliceIndex);
261 
262  if (sliceDimension > 2)
263  return;
264  if (sliceIndex >= m_SegmentationCountInSlice[timeStep][sliceDimension].size())
265  return;
266 
267  unsigned int dim0(options.dim0);
268  unsigned int dim1(options.dim1);
269 
270  int numberOfPixels(0); // number of pixels in this slice that are not 0
271 
272  unsigned int dim0max = m_SegmentationCountInSlice[timeStep][dim0].size();
273  unsigned int dim1max = m_SegmentationCountInSlice[timeStep][dim1].size();
274 
275  // scan the slice from two directions
276  // and set the flags for the two dimensions of the slice
277  for (unsigned int v = 0; v < dim1max; ++v)
278  {
279  for (unsigned int u = 0; u < dim0max; ++u)
280  {
281  DATATYPE value = *(pixelData + u + v * dim0max);
282 
283  assert((signed)m_SegmentationCountInSlice[timeStep][dim0][u] + (signed)value >=
284  0); // just for debugging. This must always be true, otherwise some counting is going wrong
285  assert((signed)m_SegmentationCountInSlice[timeStep][dim1][v] + (signed)value >= 0);
286 
287  m_SegmentationCountInSlice[timeStep][dim0][u] =
288  static_cast<unsigned int>(m_SegmentationCountInSlice[timeStep][dim0][u] + value);
289  m_SegmentationCountInSlice[timeStep][dim1][v] =
290  static_cast<unsigned int>(m_SegmentationCountInSlice[timeStep][dim1][v] + value);
291  numberOfPixels += static_cast<int>(value);
292  }
293  }
294 
295  // flag for the dimension of the slice itself
296  assert((signed)m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] + numberOfPixels >= 0);
297  m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] += numberOfPixels;
298 
299  // MITK_INFO << "scan t=" << timeStep << " from (0,0) to (" << dim0max << "," << dim1max << ") (" << pixelData << "-"
300  // << pixelData+dim0max*dim1max-1 << ") in slice " << sliceIndex << " found " << numberOfPixels << " pixels" <<
301  // std::endl;
302 }
303 
304 template <typename TPixel, unsigned int VImageDimension>
305 void mitk::SegmentationInterpolationController::ScanChangedVolume(const itk::Image<TPixel, VImageDimension> *diffImage,
306  unsigned int timeStep)
307 {
308  typedef itk::ImageSliceConstIteratorWithIndex<itk::Image<TPixel, VImageDimension>> IteratorType;
309 
310  IteratorType iter(diffImage, diffImage->GetLargestPossibleRegion());
311  iter.SetFirstDirection(0);
312  iter.SetSecondDirection(1);
313 
314  int numberOfPixels(0); // number of pixels in this slice that are not 0
315 
316  typename IteratorType::IndexType index;
317  unsigned int x = 0;
318  unsigned int y = 0;
319  unsigned int z = 0;
320 
321  iter.GoToBegin();
322  while (!iter.IsAtEnd())
323  {
324  while (!iter.IsAtEndOfSlice())
325  {
326  while (!iter.IsAtEndOfLine())
327  {
328  index = iter.GetIndex();
329 
330  x = index[0];
331  y = index[1];
332  z = index[2];
333 
334  TPixel value = iter.Get();
335 
336  assert((signed)m_SegmentationCountInSlice[timeStep][0][x] + (signed)value >=
337  0); // just for debugging. This must always be true, otherwise some counting is going wrong
338  assert((signed)m_SegmentationCountInSlice[timeStep][1][y] + (signed)value >= 0);
339 
340  m_SegmentationCountInSlice[timeStep][0][x] =
341  static_cast<unsigned int>(m_SegmentationCountInSlice[timeStep][0][x] + value);
342  m_SegmentationCountInSlice[timeStep][1][y] =
343  static_cast<unsigned int>(m_SegmentationCountInSlice[timeStep][1][y] + value);
344 
345  numberOfPixels += static_cast<int>(value);
346 
347  ++iter;
348  }
349  iter.NextLine();
350  }
351  assert((signed)m_SegmentationCountInSlice[timeStep][2][z] + numberOfPixels >= 0);
352  m_SegmentationCountInSlice[timeStep][2][z] += numberOfPixels;
353  numberOfPixels = 0;
354 
355  iter.NextSlice();
356  }
357 }
358 
359 template <typename DATATYPE>
361  const Image *volume,
362  unsigned int timeStep)
363 {
364  if (!volume)
365  return;
366  if (timeStep >= m_SegmentationCountInSlice.size())
367  return;
368 
369  ImageReadAccessor readAccess(volume, volume->GetVolumeData(timeStep));
370 
371  for (unsigned int slice = 0; slice < volume->GetDimension(2); ++slice)
372  {
373  const DATATYPE *rawVolume =
374  static_cast<const DATATYPE *>(readAccess.GetData()); // we again promise not to change anything, we'll just count
375  const DATATYPE *rawSlice = rawVolume + (volume->GetDimension(0) * volume->GetDimension(1) * slice);
376 
377  ScanChangedSlice<DATATYPE>(nullptr, SetChangedSliceOptions(2, slice, 0, 1, timeStep, rawSlice));
378  }
379 }
380 
382 {
383  unsigned int timeStep(0); // if needed, put a loop over time steps around everyting, but beware, output will be long
384 
385  MITK_INFO << "Interpolator status (timestep 0): dimensions " << m_SegmentationCountInSlice[timeStep][0].size() << " "
386  << m_SegmentationCountInSlice[timeStep][1].size() << " " << m_SegmentationCountInSlice[timeStep][2].size()
387  << std::endl;
388 
389  MITK_INFO << "Slice 0: " << m_SegmentationCountInSlice[timeStep][2][0] << std::endl;
390 
391  // row "x"
392  for (unsigned int index = 0; index < m_SegmentationCountInSlice[timeStep][0].size(); ++index)
393  {
394  if (m_SegmentationCountInSlice[timeStep][0][index] > 0)
395  MITK_INFO << "O";
396  else
397  MITK_INFO << ".";
398  }
399  MITK_INFO << std::endl;
400 
401  // rows "y" and "z" (diagonal)
402  for (unsigned int index = 1; index < m_SegmentationCountInSlice[timeStep][1].size(); ++index)
403  {
404  if (m_SegmentationCountInSlice[timeStep][1][index] > 0)
405  MITK_INFO << "O";
406  else
407  MITK_INFO << ".";
408 
409  if (m_SegmentationCountInSlice[timeStep][2].size() > index) // if we also have a z value here, then print it, too
410  {
411  for (unsigned int indent = 1; indent < index; ++indent)
412  MITK_INFO << " ";
413 
414  if (m_SegmentationCountInSlice[timeStep][2][index] > 0)
415  MITK_INFO << m_SegmentationCountInSlice[timeStep][2][index]; //"O";
416  else
417  MITK_INFO << ".";
418  }
419 
420  MITK_INFO << std::endl;
421  }
422 
423  // z indices that are larger than the biggest y index
424  for (unsigned int index = m_SegmentationCountInSlice[timeStep][1].size();
425  index < m_SegmentationCountInSlice[timeStep][2].size();
426  ++index)
427  {
428  for (unsigned int indent = 0; indent < index; ++indent)
429  MITK_INFO << " ";
430 
431  if (m_SegmentationCountInSlice[timeStep][2][index] > 0)
432  MITK_INFO << m_SegmentationCountInSlice[timeStep][2][index]; //"O";
433  else
434  MITK_INFO << ".";
435 
436  MITK_INFO << std::endl;
437  }
438 }
439 
441  unsigned int sliceIndex,
442  const mitk::PlaneGeometry *currentPlane,
443  unsigned int timeStep)
444 {
445  if (m_Segmentation.IsNull())
446  return nullptr;
447 
448  if (!currentPlane)
449  {
450  return nullptr;
451  }
452 
453  if (timeStep >= m_SegmentationCountInSlice.size())
454  return nullptr;
455  if (sliceDimension > 2)
456  return nullptr;
457  unsigned int upperLimit = m_SegmentationCountInSlice[timeStep][sliceDimension].size();
458  if (sliceIndex >= upperLimit - 1)
459  return nullptr; // can't interpolate first and last slice
460  if (sliceIndex < 1)
461  return nullptr;
462 
463  if (m_SegmentationCountInSlice[timeStep][sliceDimension][sliceIndex] > 0)
464  return nullptr; // slice contains a segmentation, won't interpolate anything then
465 
466  unsigned int lowerBound(0);
467  unsigned int upperBound(0);
468  bool bounds(false);
469 
470  for (lowerBound = sliceIndex - 1; /*lowerBound >= 0*/; --lowerBound)
471  {
472  if (m_SegmentationCountInSlice[timeStep][sliceDimension][lowerBound] > 0)
473  {
474  bounds = true;
475  break;
476  }
477 
478  if (lowerBound == 0)
479  break; // otherwise overflow and start at something like 4294967295
480  }
481 
482  if (!bounds)
483  return nullptr;
484 
485  bounds = false;
486  for (upperBound = sliceIndex + 1; upperBound < upperLimit; ++upperBound)
487  {
488  if (m_SegmentationCountInSlice[timeStep][sliceDimension][upperBound] > 0)
489  {
490  bounds = true;
491  break;
492  }
493  }
494 
495  if (!bounds)
496  return nullptr;
497 
498  // ok, we have found two neighboring slices with segmentations (and we made sure that the current slice does NOT
499  // contain anything
500  // MITK_INFO << "Interpolate in timestep " << timeStep << ", dimension " << sliceDimension << ": estimate slice " <<
501  // sliceIndex << " from slices " << lowerBound << " and " << upperBound << std::endl;
502 
503  mitk::Image::Pointer lowerMITKSlice;
504  mitk::Image::Pointer upperMITKSlice;
505  mitk::Image::Pointer resultImage;
506 
507  try
508  {
509  // Setting up the ExtractSliceFilter
511  extractor->SetInput(m_Segmentation);
512  extractor->SetTimeStep(timeStep);
513  extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep));
514  extractor->SetVtkOutputRequest(false);
515 
516  // Reslicing the current plane
517  extractor->SetWorldGeometry(currentPlane);
518  extractor->Modified();
519  extractor->Update();
520  resultImage = extractor->GetOutput();
521  resultImage->DisconnectPipeline();
522 
523  // Creating PlaneGeometry for lower slice
524  mitk::PlaneGeometry::Pointer reslicePlane = currentPlane->Clone();
525 
526  // Transforming the current origin so that it matches the lower slice
527  mitk::Point3D origin = currentPlane->GetOrigin();
528  m_Segmentation->GetSlicedGeometry(timeStep)->WorldToIndex(origin, origin);
529  origin[sliceDimension] = lowerBound;
530  m_Segmentation->GetSlicedGeometry(timeStep)->IndexToWorld(origin, origin);
531  reslicePlane->SetOrigin(origin);
532 
533  // Extract the lower slice
534  extractor = ExtractSliceFilter::New();
535  extractor->SetInput(m_Segmentation);
536  extractor->SetTimeStep(timeStep);
537  extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep));
538  extractor->SetVtkOutputRequest(false);
539 
540  extractor->SetWorldGeometry(reslicePlane);
541  extractor->Modified();
542  extractor->Update();
543  lowerMITKSlice = extractor->GetOutput();
544  lowerMITKSlice->DisconnectPipeline();
545 
546  // Transforming the current origin so that it matches the upper slice
547  m_Segmentation->GetSlicedGeometry(timeStep)->WorldToIndex(origin, origin);
548  origin[sliceDimension] = upperBound;
549  m_Segmentation->GetSlicedGeometry(timeStep)->IndexToWorld(origin, origin);
550  reslicePlane->SetOrigin(origin);
551 
552  // Extract the upper slice
553  extractor = ExtractSliceFilter::New();
554  extractor->SetInput(m_Segmentation);
555  extractor->SetTimeStep(timeStep);
556  extractor->SetResliceTransformByGeometry(m_Segmentation->GetTimeGeometry()->GetGeometryForTimeStep(timeStep));
557  extractor->SetVtkOutputRequest(false);
558 
559  extractor->SetWorldGeometry(reslicePlane);
560  extractor->Modified();
561  extractor->Update();
562  upperMITKSlice = extractor->GetOutput();
563  upperMITKSlice->DisconnectPipeline();
564 
565  if (lowerMITKSlice.IsNull() || upperMITKSlice.IsNull())
566  return nullptr;
567  }
568  catch (const std::exception &e)
569  {
570  MITK_ERROR << "Error in 2D interpolation: " << e.what();
571  return nullptr;
572  }
573 
574  // interpolation algorithm gets some inputs
575  // two segmentations (guaranteed to be of the same data type, but no special data type guaranteed)
576  // orientation (sliceDimension) of the segmentations
577  // position of the two slices (sliceIndices)
578  // one volume image (original patient image)
579  //
580  // interpolation algorithm can use e.g. itk::ImageSliceConstIteratorWithIndex to
581  // inspect the original patient image at appropriate positions
582 
585  return algorithm->Interpolate(lowerMITKSlice.GetPointer(),
586  lowerBound,
587  upperMITKSlice.GetPointer(),
588  upperBound,
589  sliceIndex,
590  sliceDimension,
591  resultImage,
592  timeStep,
593  m_ReferenceImage);
594 }
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
itk::SmartPointer< Self > Pointer
#define MITK_INFO
Definition: mitkLogMacros.h:22
#define MITK_ERROR
Definition: mitkLogMacros.h:24
void SetChangedSlice(const Image *sliceDiff, unsigned int sliceDimension, unsigned int sliceIndex, unsigned int timeStep)
Update after changing a single slice.
const void * GetData() const
Gives const access to the data.
static Pointer New()
virtual ImageDataItemPointer GetVolumeData(int t=0, int n=0, void *data=nullptr, ImportMemoryManagementType importMemoryManagement=CopyMemory) const
Definition: mitkImage.cpp:330
Pointer Clone() const
void BlockModified(bool)
Block reaction to an images Modified() events.
static SegmentationInterpolationController * GetInstance()
Get existing instance or create a new one.
void ScanChangedVolume(const itk::Image< TPixel, VImageDimension > *, unsigned int timeStep)
#define MITK_WARN
Definition: mitkLogMacros.h:23
Protected class of mitk::SegmentationInterpolationController. Don't use (you shouldn't be able to do ...
void SetChangedVolume(const Image *sliceDiff, unsigned int timeStep)
#define AccessFixedDimensionByItk_2(mitkImage, itkImageTypeFunction, dimension, arg1, arg2)
Image class for storing images.
Definition: mitkImage.h:76
void SetSegmentationVolume(const Image *segmentation)
Initialize with a whole volume.
static SegmentationInterpolationController * InterpolatorForImage(const Image *)
Find interpolator for a given image.
Image::Pointer Interpolate(unsigned int sliceDimension, unsigned int sliceIndex, const mitk::PlaneGeometry *currentPlane, unsigned int timeStep)
Generates an interpolated image for the given slice.
void ScanWholeVolume(const itk::Image< DATATYPE, 3 > *, const Image *volume, unsigned int timeStep)
std::map< const Image *, SegmentationInterpolationController * > InterpolatorMapType
#define AccessFixedDimensionByItk_1(mitkImage, itkImageTypeFunction, dimension, arg1)
void SetReferenceVolume(const Image *segmentation)
Set a reference image (original patient image) - optional.
void ScanChangedSlice(const itk::Image< DATATYPE, 2 > *, const SetChangedSliceOptions &options)
internal scan of a single slice
Describes a two-dimensional, rectangular plane.
unsigned int GetDimension() const
Get dimension of the image.
Definition: mitkImage.cpp:110
ImageReadAccessor class to get locked read access for a particular image part.
static Pointer New()
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.