Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkRandomParcellationGenerator.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 
17 //#include <mitkRandomParcellationGenerator.h>
18 
19 //To sort
20 #include <algorithm>
21 
22 //To use sqrt
23 #include <cmath>
24 
25 //Itk Iterators
26 #include <itkNeighborhoodIterator.h>
27 #include <itkImageRegionIterator.h>
28 
29 //for the use of abs()
30 #include <cstdlib>
31 
32 //Better Random Function xrand (extended range)
33 #define XRAND_MAX (RAND_MAX*(RAND_MAX + 2))
34 unsigned int xrand (void)
35 {
36  return rand () * (RAND_MAX + 1) + rand ();
37 }
38 
39 //Definition of the Set-Functions
40 template <typename TPixel, unsigned int VImageDimension>
41 void mitk::RandomParcellationGenerator<TPixel, VImageDimension>::SetImage(itk::Image<TPixel, VImageDimension> * image)
42 {
43  m_Image = image;
44 }
45 
46 template <typename TPixel, unsigned int VImageDimension>
48 {
49  m_NumberNodes = numberNodes;
50 }
51 
52 template <typename TPixel, unsigned int VImageDimension>
53 void mitk::RandomParcellationGenerator<TPixel, VImageDimension>::SetVariablesForMerging (int givenSizeOfSmallestRegion, int desiredNumberOfParcels, int givenSizeOfSmallestRegionBeginning)
54 {
55  m_GivenSizeOfSmallestRegion = givenSizeOfSmallestRegion;
56  m_DesiredNumberOfParcels = desiredNumberOfParcels;
57  m_GivenSizeOfSmallestRegionBeginning = givenSizeOfSmallestRegionBeginning;
58 }
59 
60 template <typename TPixel, unsigned int VImageDimension>
61 void mitk::RandomParcellationGenerator<TPixel, VImageDimension>::SetBoolsForMerging(bool mergingWithNumberParcels, bool mergingWithSmallestParcel, bool justMergeSmallParcels)
62 {
63  m_MergingWithNumberParcels = mergingWithNumberParcels;
64  m_MergingWithSmallestParcel = mergingWithSmallestParcel;
65  m_JustMergeSmallParcels = justMergeSmallParcels;
66 }
67 
68 //Definition of other functions
69 
70 //Calculates the Center of mass (COM) and m_SizeOfRegions
71 template <typename TPixel, unsigned int VImageDimension>
72 std::vector<double> mitk::RandomParcellationGenerator<TPixel, VImageDimension>::GetCenterOfMass( itk::ImageRegionIterator<ImageType> it_region, int valueOfRegion, bool getSizeOfRegions)
73 {
74  //Counts all tagged voxels in this region
75 
76  typedef itk::Image< TPixel, VImageDimension > ImageType;
77  int currentSizeOfRegion (0);
78  std::vector<typename ImageType::IndexType> indexVoxel;
79  std::vector<double> centerOfMass;
80  double xValue(0);
81  double yValue(0);
82  double zValue(0);
83 
84  for (it_region.GoToBegin(); !it_region.IsAtEnd(); ++it_region)
85  {
86  if(it_region.Value() == valueOfRegion)
87  {
88  indexVoxel.push_back(it_region.GetIndex());
89  currentSizeOfRegion++;
90  }
91  }
92 
93  if (getSizeOfRegions == true)
94  {
95  m_SizeOfRegions.push_back(currentSizeOfRegion);
96  }
97 
98  xValue = 0;
99  yValue = 0;
100  zValue = 0;
101 
102  for (int i = 0; i < currentSizeOfRegion; i++)
103  {
104  xValue += indexVoxel[i][0];
105  yValue += indexVoxel[i][1];
106  zValue += indexVoxel[i][2];
107  }
108 
109  centerOfMass.push_back(xValue/currentSizeOfRegion);
110  centerOfMass.push_back(yValue/currentSizeOfRegion);
111  centerOfMass.push_back(zValue/currentSizeOfRegion);
112 
113  return centerOfMass;
114 }
115 
116 template <typename TPixel, unsigned int VImageDimension>
117 double mitk::RandomParcellationGenerator<TPixel, VImageDimension>::GetDistance( std::vector<double> centerOfMass, typename itk::Image< TPixel, VImageDimension >::IndexType indexNewVoxel)
118 {
119  double distancetemp(0);
120 
121  for (int i = 0; i < 3; i++)
122  {
123  distancetemp += (centerOfMass[i] - indexNewVoxel[i])*(centerOfMass[i] - indexNewVoxel[i]);
124  }
125  return sqrt(distancetemp);
126 }
127 
128 template <typename TPixel, unsigned int VImageDimension>
129 double mitk::RandomParcellationGenerator<TPixel, VImageDimension>::GetDistanceVector( std::vector<double> centerOfMass, std::vector<double> indexNewVoxel)
130 {
131  double distancetemp(0);
132 
133  for (int i = 0; i < 3; i++)
134  {
135  distancetemp += (centerOfMass[i] - indexNewVoxel[i])*(centerOfMass[i] - indexNewVoxel[i]);
136  }
137  return sqrt(distancetemp);
138 }
139 
140 
141 template < typename TPixel, unsigned int VImageDimension >
143 {
144  int distancemin;
145  distancemin = distance[0];
146 
147  for (int i = 1; i < distance.size(); i++)
148  {
149  if (distance[i] < distancemin)
150  {
151  distancemin = distance[i];
152  }
153  }
154  return distancemin;
155 }
156 
157 template < typename TPixel, unsigned int VImageDimension >
159 {
160  double distancemin;
161  distancemin = distance[0];
162 
163  for (int i = 1; i < distance.size(); i++)
164  {
165  if (distance[i] < distancemin)
166  {
167  distancemin = distance[i];
168  }
169  }
170  return distancemin;
171 }
172 
173 template < typename TPixel, unsigned int VImageDimension >
175 {
176  typedef itk::Image< TPixel, VImageDimension > TemplateImageType;
177 
178  itk::ImageRegionIterator<TemplateImageType> it_image(m_Image, m_Image->GetLargestPossibleRegion());
179 
180  TPixel minimumValue, maximumValue;
181 
182  it_image.GoToBegin();
183  maximumValue = minimumValue = it_image.Value();
184 
185  for(it_image.GoToBegin(); !it_image.IsAtEnd(); ++it_image)
186  {
187  if ( it_image.Value() < minimumValue )
188  {
189  minimumValue = it_image.Value();
190  }
191  else
192  {
193  if ( it_image.Value() > maximumValue )
194  {
195  maximumValue = it_image.Value();
196  }
197  }
198  }
199 
200  int range = int ( maximumValue - minimumValue ); //needs to be castable to int
201  int offset = int ( minimumValue );
202 
203  if ( range < 0 ) //error
204  {
205  return;
206  }
207 
208  std::vector< unsigned int > histogram;
209  histogram.resize( range + 1, 0 );
210 
211  for(it_image.GoToBegin(); !it_image.IsAtEnd(); ++it_image)
212  {
213  histogram[ int ( it_image.Value() ) - offset ] += 1;
214  }
215 
216  int gapCounter = 0; //this variable will be used to count the empty labels
217 
218  //stores how much has to be subtracted from the image to remove gaps
219  std::vector< TPixel > subtractionStorage;
220  subtractionStorage.resize( range + 1, 0 );
221 
222  for( int index = 0; index <= range; index++ )
223  {
224  if( histogram[ index ] == 0 )
225  {
226  gapCounter++; //if the label is empty, increase gapCounter
227  }
228  else
229  {
230  subtractionStorage[ index ] = TPixel ( gapCounter );
231  }
232  }
233 
234  //remove gaps from label image
235  for(it_image.GoToBegin(); !it_image.IsAtEnd(); ++it_image)
236  {
237  it_image.Value() = it_image.Value() - subtractionStorage[ int ( it_image.Value() ) ];
238  }
239 }
240 
241 
242 
243 template < typename TPixel, unsigned int VImageDimension >
244 typename itk::Image< TPixel, VImageDimension >::RegionType mitk::RandomParcellationGenerator<TPixel, VImageDimension>::ExtendedRegion(typename itk::Image< TPixel, VImageDimension >::RegionType chosenRegion, typename ImageType::IndexType indexChosenVoxel)
245 {
246  itk::Size<3> size = chosenRegion.GetSize();
247  itk::Index<3> start = chosenRegion.GetIndex();
248  std::vector<int> greatestValues;
249  std::vector<int> smallestValues;
250 
251  for(int i = 0; i < 3; i++)
252  {
253  greatestValues.push_back(start[i] + size[i]);
254  smallestValues.push_back(start[i]);
255  }
256 
257  for(int i = 0; i < 3; i++)
258  {
259  if (indexChosenVoxel[i] == greatestValues[i])
260  {
261  size[i] += 1;
262  }
263  if (indexChosenVoxel[i] == smallestValues[i] - 1)
264  {
265  start[i] -= 1;
266  size[i] += 1;
267  }
268  }
269 
270  chosenRegion.SetSize(size);
271  chosenRegion.SetIndex(start);
272 
273  return chosenRegion;
274 }
275 
276 template < typename TPixel, unsigned int VImageDimension >
277 typename itk::Image< TPixel, VImageDimension >::RegionType mitk::RandomParcellationGenerator<TPixel, VImageDimension>::ExtendedRegionNotNeighbor(typename itk::Image< TPixel, VImageDimension >::RegionType chosenRegion, typename itk::Image< TPixel, VImageDimension >::RegionType smallestRegion)
278 {
279  itk::Size<3> sizeChosenRegion = chosenRegion.GetSize();
280  itk::Index<3> startChosenRegion = chosenRegion.GetIndex();
281  itk::Size<3> sizeSmallestRegion = smallestRegion.GetSize();
282  itk::Index<3> startSmallestRegion = smallestRegion.GetIndex();
283 
284  //size and start of the new merged region
285  itk::Size<3> size;
286  itk::Index<3> start;
287 
288  //Calculate the size of the merged region
289  for (int i = 0; i < 3; i++)
290  {
291  //Overlapping (start of ChosenRegion before start of SmallestRegion)
292  if ((startChosenRegion[i] + sizeChosenRegion[i] > startSmallestRegion[i]) &&
293  (startChosenRegion[i] <= startSmallestRegion[i]))
294  {
295  //Not integrated
296  if (startSmallestRegion[i] + sizeSmallestRegion[i] > startChosenRegion[i] + sizeChosenRegion[i])
297  {
298  size[i] = abs(startSmallestRegion[i] - startChosenRegion[i]) + sizeSmallestRegion[i];
299  }
300  //Integrated
301  else
302  {
303  size[i] = std::max(sizeChosenRegion[i], sizeSmallestRegion[i]);
304  }
305  }
306 
307  //Overlapping (start of SmallestRegion before start of ChosenRegion)
308  if ((startSmallestRegion[i] + sizeSmallestRegion[i] > startChosenRegion[i]) &&
309  (startChosenRegion[i] >= startSmallestRegion[i]))
310  {
311  //Not integrated
312  if (startChosenRegion[i] + sizeChosenRegion[i] > startSmallestRegion[i] + sizeSmallestRegion[i])
313  {
314  size[i] = abs(startSmallestRegion[i] - startChosenRegion[i]) + sizeChosenRegion[i];
315  }
316  //Integrated
317  else
318  {
319  size[i] = std::max(sizeChosenRegion[i], sizeSmallestRegion[i]);
320  }
321  }
322 
323  //No overlapping (right next to each other)
324  if ((startChosenRegion[i] + sizeChosenRegion[i] == startSmallestRegion[i]) ||
325  (startSmallestRegion[i] + sizeSmallestRegion[i] == startChosenRegion[i]))
326  {
327  size[i] = sizeChosenRegion[i] + sizeSmallestRegion[i];
328  }
329 
330  //Isolated
331  if ((startChosenRegion[i] + sizeChosenRegion[i] < startSmallestRegion[i]) ||
332  (startSmallestRegion[i] + sizeSmallestRegion[i] < startChosenRegion[i]))
333  {
334  if(startChosenRegion[i] + sizeChosenRegion[i] < startSmallestRegion[i])
335  {
336  size[i] = abs(startChosenRegion[i] - startSmallestRegion[i]) +1;
337  }
338  if(startSmallestRegion[i] + sizeSmallestRegion[i] < startChosenRegion[i])
339  {
340  size[i] = abs(startChosenRegion[i] - startSmallestRegion[i]) + sizeChosenRegion[i];
341  }
342  }
343  }
344 
345  //Calculate the start point of the new region, choose the smallest value respectively
346  for (int i = 0; i < 3; i++)
347  {
348  start[i] = std::min(startChosenRegion[i], startSmallestRegion[i]);
349  }
350 
351  chosenRegion.SetSize(size);
352  chosenRegion.SetIndex(start);
353 
354  return chosenRegion;
355 }
356 
357 //checks if a number is already an element of the vector
358 template < typename TPixel, unsigned int VImageDimension >
360 {
361  for (int i = 0; i < vec.size(); i++)
362  {
363  if (vec[i] == number)
364  {
365  return false;
366  }
367  }
368  return true;
369 }
370 
371 template < typename TPixel, unsigned int VImageDimension >
373 {
374  typedef itk::Image< TPixel, VImageDimension > ImageType;
375 
377  counter.SetImage(m_Image);
378  counter.SetRegion(m_Image->GetLargestPossibleRegion());
379 
380  int numberVoxels = counter.VoxelWithValue(1);
381 
382  //Create vector with unique numbers
383 
384  std::vector<int> randomvector;
385  int randomnumber;
386 
387  //for-loop to get (m_NumberNodes) randomly chosen seed-points
388  for (int j = 0; j < m_NumberNodes; j++)
389  {
390  //xrand() is the same as rand() but with an extended range
391  randomnumber = xrand() % numberVoxels +1;
392 
393  //Bool-function (IsUnique)
394  if (this->IsUnique(randomnumber, randomvector))
395  {
396  randomvector.push_back(randomnumber);
397  }
398  else
399  {
400  j--;
401  }
402  }
403 
404  //sorts the numbers in randomvector according to their size
405  sort (randomvector.begin(), randomvector.end());
406 
407  //assign new values to the chosen voxels and build the regions
408 
409  itk::Size<3> originalSize;
410  originalSize.Fill(1);
411  itk::Index<3> originalStart;
412  typename ImageType::IndexType currentIndex;
413 
414  int valueCounter(0);
415  int regionNumber(0);
416  typename ImageType::RegionType tempRegion;
417  itk::ImageRegionIterator<ImageType> it_image(m_Image, m_Image->GetLargestPossibleRegion());
418 
419 
420  for (it_image.GoToBegin(); !it_image.IsAtEnd(), regionNumber < m_NumberNodes; ++it_image)
421  {
422  if (it_image.Value() >= 1)
423  {
424  valueCounter++;
425  if (valueCounter == randomvector[regionNumber])
426  {
427  it_image.Value() = regionNumber+3;
428 
429  //Build the Regions
430  currentIndex = it_image.GetIndex();
431  originalStart[0] = currentIndex[0];
432  originalStart[1] = currentIndex[1];
433  originalStart[2] = currentIndex[2];
434 
435  tempRegion.SetSize(originalSize);
436  tempRegion.SetIndex(originalStart);
437  std::pair<RegionType, int> tempPair (tempRegion, regionNumber+3);
438  m_OddRegions.push_back(tempPair);
439  regionNumber++;
440  }
441  }
442  }
443 }
444 
445 template < typename TPixel, unsigned int VImageDimension >
447 {
448  typedef itk::Image< int, VImageDimension > IntegerImageType;
449  typedef itk::NeighborhoodIterator< IntegerImageType > NeighborhoodIteratorType;
450  typename NeighborhoodIteratorType::RadiusType radius;
451  radius.Fill(1);
452 
453  std::pair<RegionType, int> chosenRegion;
454  int thisRegion;
455  typename IntegerImageType::IndexType indexChosenVoxel;
456  std::vector<double> distance;
457  std::vector<int> possibleNeighbor;
458  int thisNeighbor;
459  std::vector<double> centerOfMass;
460  typename IntegerImageType::IndexType indexNewVoxel;
461  std::vector<typename IntegerImageType::IndexType> indexNewVoxelVector;
462  int numberAvailableVoxels(0);
463  bool validVoxelNotYetFound;
464  bool checkRegionsOfOddSize(true);
465  std::vector<std::pair<RegionType,int> > * possibleRegions;
466  int costValue;
467 
468 
469  //As long as there are still non invalid regions (of even or odd size) add a voxel to them in each iteration step
470  while(m_OddRegions.size() != 0 || m_EvenRegions.size() != 0)
471  {
472  //When we change from odd to even or vice versa
473  if (m_OddRegions.size() == 0)
474  {
475  checkRegionsOfOddSize = false;
476  possibleRegions = &m_EvenRegions;
477  }
478 
479  if (m_EvenRegions.size() == 0)
480  {
481  checkRegionsOfOddSize = true;
482  possibleRegions = &m_OddRegions;
483  }
484 
485  //All our smallest regions are contained in possibleRegions, choose one randomly and check whether we can expand it
486 
487  thisRegion = xrand() % possibleRegions->size();
488  chosenRegion = (*possibleRegions)[thisRegion];
489 
490  //Calculate Center of mass (COM)
491  centerOfMass.clear();
492  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, chosenRegion.first);
493  centerOfMass = this->GetCenterOfMass(it_region, chosenRegion.second, false);
494 
495  numberAvailableVoxels = 0;
496  NeighborhoodIteratorType it_neighborhood(radius, m_Image, chosenRegion.first);
497 
498  //Check if there are any valid voxels which we can add to our chosen region
499 
500  for (it_neighborhood.GoToBegin(); !it_neighborhood.IsAtEnd(); ++it_neighborhood)
501  {
502  if (it_neighborhood.GetCenterPixel() == chosenRegion.second)
503  {
504  //Now the Neighborhood of a Parcel-Voxel is available
505  //Check if there is a valid voxel
506 
507  for (int k = 0; k <27; k++)
508  {
509  if (it_neighborhood.GetPixel(k) == 1 && (k == 4 || k == 10 || k == 12 || k == 14 || k == 16 || k == 22) )
510  {
511  // Calculate the distance btw. the centerOfMass and the NewVoxel
512  indexNewVoxel = it_neighborhood.GetIndex(k);
513  indexNewVoxelVector.push_back(indexNewVoxel);
514 
515  distance.push_back(this->GetDistance(centerOfMass, indexNewVoxel));
516 
517  numberAvailableVoxels++;
518  }
519  }
520  }
521  }
522 
523  if (numberAvailableVoxels > 0)
524  {
525  validVoxelNotYetFound = true;
526  while (validVoxelNotYetFound && distance.size() > 0)
527  {
528  possibleNeighbor.clear();
529  //Find the possible neighbors with the smallest distance to the COM of the region
530 
531  for (int i = 0; i < distance.size(); i++)
532  {
533  if (distance[i] == this->SmallestValue(distance))
534  {
535  possibleNeighbor.push_back(i);
536  }
537  }
538 
539  //Choose a random voxel and check if it is valid
540  //Get the index of this voxel
541  thisNeighbor = xrand() % possibleNeighbor.size();
542  indexChosenVoxel = indexNewVoxelVector[possibleNeighbor[thisNeighbor]];
543 
544  //Check if we now have to expand the region due to the possible new voxel
545  std::pair<RegionType, int> chosenRegionChanged;
546  chosenRegionChanged.first = this->ExtendedRegion(chosenRegion.first, indexChosenVoxel);
547  chosenRegionChanged.second = chosenRegion.second;
548 
549  //Constraints
550 
552  costCalculator.SetImage(m_Image);
553  costCalculator.SetRegion(chosenRegionChanged);
554  costValue = costCalculator.CalculateCost();
555 
556 
557  if(costValue >= costCalculator.MaximalValue())//Constraints are fulfilled
558  {
559  validVoxelNotYetFound = false;
560  chosenRegion = chosenRegionChanged;
561  }
562  else //Constraints are not fulfilled
563  {
564  validVoxelNotYetFound = true; //So we are still in the while loop
565  distance.erase(distance.begin() + possibleNeighbor[thisNeighbor]);
566  indexNewVoxelVector.erase(indexNewVoxelVector.begin() + possibleNeighbor[thisNeighbor]);
567 
568  if (distance.size() == 0)
569  {
570  if (checkRegionsOfOddSize == true)
571  {
572  m_InvalidRegions.push_back(chosenRegion);
573  m_OddRegions.erase(m_OddRegions.begin() + thisRegion);
574  }
575  else
576  {
577  m_InvalidRegions.push_back(chosenRegion);
578  m_EvenRegions.erase(m_EvenRegions.begin() + thisRegion);
579  }
580  }
581  }
582  }
583 
584 
585  if (distance.size() > 0)
586  {
587  //Change the value of the new voxel
588  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_Image->GetLargestPossibleRegion());
589  it_region.SetIndex(indexChosenVoxel);
590  it_region.Value() = chosenRegion.second;
591 
592  //Change from even to odd or vice versa
593  if (checkRegionsOfOddSize == true)
594  {
595  m_EvenRegions.push_back(chosenRegion);
596  m_OddRegions.erase(m_OddRegions.begin() + thisRegion);
597  }
598  else
599  {
600  m_OddRegions.push_back(chosenRegion);
601  m_EvenRegions.erase(m_EvenRegions.begin() + thisRegion);
602  }
603 
604  distance.clear();
605  indexNewVoxelVector.clear();
606  }
607  }
608  else
609  {
610  //The region can't be extended, put it into the m_InvalidRegions vector
611  m_InvalidRegions.push_back((*possibleRegions)[thisRegion]);
612  possibleRegions->erase(possibleRegions->begin() + thisRegion);
613  }
614  }
615 }
616 
617 template < typename TPixel, unsigned int VImageDimension >
619 {
620  typedef itk::Image< int, VImageDimension > IntegerImageType;
621  typedef itk::NeighborhoodIterator< IntegerImageType > NeighborhoodIteratorType;
622  typename NeighborhoodIteratorType::RadiusType radius;
623  radius.Fill(1);
624  NeighborhoodIteratorType it_neighborhood(radius, m_Image, m_Image->GetLargestPossibleRegion());
625 
626  std::vector<int> valueOfRegions;
627  std::vector<int> sizeOfPossibleRegions;
628  std::vector<int> valueOfPossibleRegions;
629  std::vector<int> indexInvalidRegions;
630  std::vector<int> indexOfPossibleInvalidRegions;
631  std::vector<double> distance;
632  std::vector<int> possibleNeighbor;
633  typename IntegerImageType::IndexType indexNewVoxel;
634  std::vector<double> centerOfMass;
635  int thisNeighbor;
636  std::pair<RegionType, int> chosenRegion;
637  typename IntegerImageType::IndexType indexChosenVoxel;
638 
640  wholeImageCounter.SetImage(m_Image);
641  wholeImageCounter.SetRegion(m_Image->GetLargestPossibleRegion());
642  int voxelsWithValueOne = wholeImageCounter.VoxelWithValue(1);
643 
644  if (voxelsWithValueOne > 0)
645  {
646  std::stringstream message;
647  message << "Edge " << voxelsWithValueOne << endl;
648  MITK_DEBUG << message.str();
649  for (it_neighborhood.GoToBegin(); !it_neighborhood.IsAtEnd(); ++it_neighborhood)
650  {
651  if (it_neighborhood.GetCenterPixel() == 1) //Found Voxel with Value 1
652  {
653  m_SizeOfRegions.clear();
654  distance.clear();
655  valueOfRegions.clear();
656  possibleNeighbor.clear();
657  indexInvalidRegions.clear();
658  indexOfPossibleInvalidRegions.clear();
659  indexNewVoxel = it_neighborhood.GetIndex();
660 
661  //Add this voxel to a region
662  //first check if it shares an edge or vertex with (at least) one region
663  for (int k = 0; k <27; k++)
664  {
665  if (it_neighborhood.GetPixel(k) > 2)
666  {
667  for (int i = 0; i < m_InvalidRegions.size(); i++)
668  {
669  if (it_neighborhood.GetPixel(k) == m_InvalidRegions[i].second)
670  {
671  valueOfRegions.push_back(m_InvalidRegions[i].second);
672  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_InvalidRegions[i].first);
673 
674  //Gets center of mass and the m_SizeOfRegions
675  centerOfMass = this->GetCenterOfMass(it_region, m_InvalidRegions[i].second, true);
676  distance.push_back(this->GetDistance(centerOfMass, indexNewVoxel));
677  indexInvalidRegions.push_back(i);
678  break;
679  }
680  }
681  }
682  }
683  //Minimal distance
684  if (distance.size() != 0)
685  {
686  for (int i = 0; i < distance.size(); i++)
687  {
688  if (distance[i] == this->SmallestValue(distance))
689  {
690  possibleNeighbor.push_back(i);
691  }
692  }
693  distance.clear();
694 
695  //Check the Size of the Regions and get the smallest one!
696  //First get the regions with the same distance to COM
697  for (int i = 0; i < possibleNeighbor.size(); i++)
698  {
699  sizeOfPossibleRegions.push_back(m_SizeOfRegions[possibleNeighbor[i]]);
700  valueOfPossibleRegions.push_back(valueOfRegions[possibleNeighbor[i]]);
701  indexOfPossibleInvalidRegions.push_back(indexInvalidRegions[possibleNeighbor[i]]);
702  }
703  possibleNeighbor.clear();
704  indexInvalidRegions.clear();
705 
706  for (int i = 0; i < sizeOfPossibleRegions.size(); i++)
707  {
708  if (sizeOfPossibleRegions[i] == this->SmallestValue(sizeOfPossibleRegions))
709  {
710  possibleNeighbor.push_back(i);
711  }
712  }
713  sizeOfPossibleRegions.clear();
714 
715  thisNeighbor = xrand() % possibleNeighbor.size();
716 
717  it_neighborhood.SetCenterPixel(valueOfPossibleRegions[thisNeighbor]);
718  valueOfPossibleRegions.clear();
719 
720  //region extansion if necessary
721  //first we need to get the right region
722  chosenRegion = m_InvalidRegions[indexOfPossibleInvalidRegions[thisNeighbor]];
723  indexChosenVoxel = it_neighborhood.GetIndex();
724 
725  m_InvalidRegions[indexOfPossibleInvalidRegions[thisNeighbor]].first = this->ExtendedRegion(chosenRegion.first, indexChosenVoxel);
726 
727  indexOfPossibleInvalidRegions.clear();
728 
729  }
730  }
731  }
732  }
733 }
734 
735 template < typename TPixel, unsigned int VImageDimension >
737 {
738  typedef itk::Image< int, VImageDimension > IntegerImageType;
739 
740  std::vector<int> valueOfRegions;
741  std::vector<int> sizeOfPossibleRegions;
742  std::vector<int> valueOfPossibleRegions;
743  std::vector<int> indexInvalidRegions;
744  std::vector<int> indexOfPossibleInvalidRegions;
745  std::vector<double> distance;
746  std::vector<int> possibleNeighbor;
747  typename IntegerImageType::IndexType indexNewVoxel;
748  std::vector<double> centerOfMass;
749  int thisNeighbor;
750  std::pair<RegionType, int> chosenRegion;
751  typename IntegerImageType::IndexType indexChosenVoxel;
752 
754  wholeImageCounter.SetImage(m_Image);
755  wholeImageCounter.SetRegion(m_Image->GetLargestPossibleRegion());
756  int voxelsWithValueOne = wholeImageCounter.VoxelWithValue(1);
757 
758  if (voxelsWithValueOne > 0)
759  {
760  std::stringstream message;
761  message << "Isolated " << voxelsWithValueOne << endl;
762  MITK_DEBUG << message.str();
763  std::vector<std::vector<double> > comOfRegions;
764  m_SizeOfRegions.clear();
765 
766  //Calculate all center of mass
767  for (int i = 0; i < m_InvalidRegions.size(); i++)
768  {
769  valueOfRegions.push_back(m_InvalidRegions[i].second);
770  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_InvalidRegions[i].first);
771 
772  centerOfMass = this->GetCenterOfMass(it_region, m_InvalidRegions[i].second, true);
773  comOfRegions.push_back(centerOfMass);
774  }
775 
776  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_Image->GetLargestPossibleRegion());
777 
778  for (it_region.GoToBegin(); !it_region.IsAtEnd(); ++it_region)
779  {
780  if (it_region.Value() == 1) //Found Voxel with Value 1
781  {
782  indexNewVoxel = it_region.GetIndex();
783 
784  //distance calculation
785  for (int j = 0; j < m_InvalidRegions.size(); j++)
786  {
787  distance.push_back(this->GetDistance(comOfRegions[j], indexNewVoxel));
788  indexInvalidRegions.push_back(j);
789 
790  }
791 
792  //Minimal distance
793  if (distance.size() != 0)
794  {
795  for (int i = 0; i < distance.size(); i++)
796  {
797  if (distance[i] == this->SmallestValue(distance))
798  {
799  possibleNeighbor.push_back(i);
800  }
801  }
802  distance.clear();
803 
804  //Check the Size of the Regions and get the smallest one!
805  //First get the regions with the same distance to COM
806  for (int i = 0; i < possibleNeighbor.size(); i++)
807  {
808  sizeOfPossibleRegions.push_back(m_SizeOfRegions[possibleNeighbor[i]]);
809  valueOfPossibleRegions.push_back(valueOfRegions[possibleNeighbor[i]]);
810  indexOfPossibleInvalidRegions.push_back(indexInvalidRegions[possibleNeighbor[i]]);
811  }
812  possibleNeighbor.clear();
813  indexInvalidRegions.clear();
814 
815  for (int i = 0; i < sizeOfPossibleRegions.size(); i++)
816  {
817  if (sizeOfPossibleRegions[i] == this->SmallestValue(sizeOfPossibleRegions))
818  {
819  possibleNeighbor.push_back(i);
820  }
821  }
822  sizeOfPossibleRegions.clear();
823 
824  thisNeighbor = xrand() % possibleNeighbor.size();
825 
826  it_region.Value() = valueOfPossibleRegions[thisNeighbor];
827  valueOfPossibleRegions.clear();
828  possibleNeighbor.clear();
829 
830  //region extansion if necessary
831  //first we need to get the right region
832  chosenRegion = m_InvalidRegions[indexOfPossibleInvalidRegions[thisNeighbor]];
833  indexChosenVoxel = it_region.GetIndex();
834  RegionType voxelRegion;
835  voxelRegion.SetIndex(indexChosenVoxel);
836  itk::Size<3> voxelSize;
837  voxelSize.Fill(1);
838  voxelRegion.SetSize(voxelSize);
839 
840  m_InvalidRegions[indexOfPossibleInvalidRegions[thisNeighbor]].first = this->ExtendedRegionNotNeighbor(chosenRegion.first, voxelRegion);
841 
842  indexOfPossibleInvalidRegions.clear();
843  }
844  }
845  }
846  }
847 }
848 
849 template < typename TPixel, unsigned int VImageDimension >
851 {
852  typedef itk::Image< int, VImageDimension > IntegerImageType;
853  typedef itk::NeighborhoodIterator< IntegerImageType > NeighborhoodIteratorType;
854  typename NeighborhoodIteratorType::RadiusType radius;
855  radius.Fill(1);
856 
857  //Calculate the m_SizeOfRegions
858  this->GetSizeOfRegions();
859 
860  int sizeOfSmallestRegion = this->SmallestValue(m_SizeOfRegions);
861 
862  std::vector<int> smallRegions;
863  std::vector<int> indexSmallRegions;
864  int thisIndex;
865  int thisSmallRegion;
866  bool mergingStillPossible(true);
867  std::vector<int> smallDistances;
868  std::vector<double> centerOfMassPossibleRegion;
869  std::vector<std::vector<double> > comOfRegions;
870  std::vector<int> sizeSmallRegions;
871  std::vector<int> smallestRegions;
872  bool tooManyParcels(true);
873  int currentNumberOfParcels = m_NumberNodes;
874  int sizeNewParcel;
875  std::vector<double> costFunctionValue;
876  std::vector<double> smallestCostFunctionValues;
877  std::vector<int> valueOfPossibleRegions;
878  std::vector<double> distance;
879  std::pair<RegionType, int> chosenRegion;
880  std::vector<double> centerOfMass;
881  bool smallRegionFound(false);
882  bool hasNeighbors(false);
883 
884  while (sizeOfSmallestRegion < m_GivenSizeOfSmallestRegion && mergingStillPossible && tooManyParcels)
885  {
886  smallRegions.clear();
887  sizeSmallRegions.clear();
888  smallestRegions.clear();
889 
890  //Find all small Regions
891  for (int i = 0; i < m_SizeOfRegions.size(); i++)
892  {
893  if (m_SizeOfRegions[i] < m_GivenSizeOfSmallestRegion)
894  {
895  smallRegions.push_back(i);
896  sizeSmallRegions.push_back(m_SizeOfRegions[i]);
897  }
898  }
899 
900  //Find one of the smallest regions
901  for (int i = 0; i < sizeSmallRegions.size(); i++)
902  {
903  if (sizeSmallRegions[i] == this->SmallestValue(sizeSmallRegions))
904  {
905  smallestRegions.push_back(smallRegions[i]);
906  }
907  }
908 
909  if (smallestRegions.size() > 0 && this->SmallestValue(sizeSmallRegions) < m_GivenSizeOfSmallestRegionBeginning)
910  {
911  thisSmallRegion = rand() % smallestRegions.size();
912 
913  //Choose a random small region
914  //First check if it has direct neighbors
915  NeighborhoodIteratorType it_neighborhood(radius, m_Image, m_InvalidRegions[smallestRegions[thisSmallRegion]].first);
916  valueOfPossibleRegions.clear();
917  smallRegionFound = false;
918  hasNeighbors = false;
919 
920  for (it_neighborhood.GoToBegin(); !it_neighborhood.IsAtEnd(); ++it_neighborhood)
921  {
922  for (int k = 0; k <27; k++)
923  {
924  if ((k == 4 || k == 10 || k == 12 || k == 14 || k == 16 || k == 22) &&//over faces
925  it_neighborhood.GetPixel(k) > 2 //found value which belongs to a region != the chosen region
926  && it_neighborhood.GetPixel(k) != m_InvalidRegions[smallestRegions[thisSmallRegion]].second)
927  {
928  hasNeighbors = true;
929  for (int i = 0; i < smallRegions.size(); i++)
930  {
931  if (it_neighborhood.GetPixel(k) == m_InvalidRegions[smallRegions[i]].second) //The value belongs to a small region
932  {
933  valueOfPossibleRegions.push_back(it_neighborhood.GetPixel(k));
934  smallRegionFound = true;
935  break;
936  }
937  }
938  }
939  }
940  }
941 
942  if (hasNeighbors == true)
943  {
944  if (smallRegionFound == true)//the region has direct neighbors
945  {
946  //valueOfPossibleRegions may contain some values, that are the same. Erase the duplicates and sort the vector.
947  std::sort(valueOfPossibleRegions.begin(), valueOfPossibleRegions.end());
948  valueOfPossibleRegions.erase(std::unique(valueOfPossibleRegions.begin(), valueOfPossibleRegions.end()), valueOfPossibleRegions.end());
949 
950  indexSmallRegions.clear();
951 
952  //Get the values of the possible regions!
953  for (int i = 0; i < valueOfPossibleRegions.size() ; i++)
954  {
955  for (int j = 0; j < m_InvalidRegions.size(); j++)
956  {
957  if (valueOfPossibleRegions[i] == m_InvalidRegions[j].second)
958  {
959  indexSmallRegions.push_back(j);
960  break;
961  }
962  }
963  }
964 
965  //take the region with the greatest value of the cost function
966  //the cost function depends on the distance (btw the 2 COMs) and the size of the merged parcel
967 
968  //First get all the distances
969  itk::ImageRegionIterator<IntegerImageType> it_regionSmallest(m_Image, m_InvalidRegions[smallestRegions[thisSmallRegion]].first);
970  std::vector<double> centerOfMassCurrentRegion = this->GetCenterOfMass(it_regionSmallest, m_InvalidRegions[smallestRegions[thisSmallRegion]].second, false);
971 
972  distance.clear();
973 
974  for (int i = 0; i < indexSmallRegions.size(); i++)
975  {
976  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_InvalidRegions[indexSmallRegions[i]].first);
977  centerOfMassPossibleRegion = this->GetCenterOfMass(it_region, m_InvalidRegions[indexSmallRegions[i]].second, false);
978  distance.push_back(this->GetDistanceVector(centerOfMassCurrentRegion, centerOfMassPossibleRegion));
979  }
980  smallDistances.clear();
981 
982  if (m_MergingWithSmallestParcel && m_JustMergeSmallParcels)
983  {
984  //If there are small Regions with equal distance btw the COMs choose a random one
985  for (int i = 0; i < distance.size(); i++)
986  {
987  if (distance[i] == this->SmallestValue(distance))
988  {
989  smallDistances.push_back(indexSmallRegions[i]);
990  }
991  }
992  thisIndex = smallDistances[rand() % smallDistances.size()];
993  }
994 
995  //Calculate the cost function values
996 
997  if (m_MergingWithNumberParcels || (m_MergingWithSmallestParcel && !m_JustMergeSmallParcels))
998  {
999  costFunctionValue.clear();
1000 
1001  for (int i = 0; i < indexSmallRegions.size(); i++)
1002  {
1003  if (distance[i] == 0)
1004  {
1005  distance[i] = 0.01;
1006  }
1007  sizeNewParcel = (m_SizeOfRegions[smallestRegions[thisSmallRegion]] + m_SizeOfRegions[indexSmallRegions[i]]);
1008 
1009  costFunctionValue.push_back(-(1.0/distance[i]) - (1.0/sizeNewParcel)); //1, 2 or 3 at sizeNewParcel
1010  }
1011  smallestCostFunctionValues.clear();
1012 
1013  for (int i = 0; i < indexSmallRegions.size(); i++)
1014  {
1015  if (costFunctionValue[i] == this->SmallestValue(costFunctionValue))
1016  {
1017  smallestCostFunctionValues.push_back(indexSmallRegions[i]);
1018  }
1019  }
1020  thisIndex = smallestCostFunctionValues[rand() % smallestCostFunctionValues.size()];
1021  }
1022 
1023  //m_InvalidRegions[thisIndex].first is the region we want to merge with our current one
1024 
1025  //Colour the small parcel like the found one
1026  for (it_regionSmallest.GoToBegin(); !it_regionSmallest.IsAtEnd(); ++it_regionSmallest)
1027  {
1028  if (it_regionSmallest.Value() == m_InvalidRegions[smallestRegions[thisSmallRegion]].second)
1029  {
1030  it_regionSmallest.Value() = m_InvalidRegions[thisIndex].second;
1031  }
1032  }
1033 
1034  //Expand the region of the new parcel if necessary!
1035 
1036  chosenRegion = m_InvalidRegions[thisIndex];
1037  chosenRegion.first = this->ExtendedRegionNotNeighbor(chosenRegion.first, m_InvalidRegions[smallestRegions[thisSmallRegion]].first);
1038  m_InvalidRegions[thisIndex] = chosenRegion;
1039 
1040  //Erase the smallest Region from m_InvalidRegions
1041  m_InvalidRegions.erase(m_InvalidRegions.begin() + smallestRegions[thisSmallRegion]);
1042  currentNumberOfParcels--;
1043 
1044  if (currentNumberOfParcels <= m_DesiredNumberOfParcels)
1045  {
1046  tooManyParcels = false;
1047  }
1048 
1049  //m_SizeOfRegions changed, get the new one
1050  this->GetSizeOfRegions();
1051  }
1052 
1053  else //No merging possible, erase this small region from m_InvalidRegions and from m_SizeOfRegions
1054  {
1055  m_InvalidRegions.erase(m_InvalidRegions.begin() + smallestRegions[thisSmallRegion]);
1056  currentNumberOfParcels--;
1057  m_SizeOfFinishedRegions.push_back(m_SizeOfRegions[smallestRegions[thisSmallRegion]]);
1058  m_SizeOfRegions.erase(m_SizeOfRegions.begin() + smallestRegions[thisSmallRegion]);
1059 
1060  if (currentNumberOfParcels <= m_DesiredNumberOfParcels)
1061  {
1062  tooManyParcels = false;
1063  }
1064  }
1065  }
1066  else //There are no 6-connected neighborhood-regions. Try to merge this small region to a near region (according to the Center of Mass)
1067  {
1068  comOfRegions.clear();
1069  //Calculate the center of mass of all small regions
1070  for (int i = 0; i < smallRegions.size(); i++)
1071  {
1072  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_InvalidRegions[smallRegions[i]].first);
1073 
1074  centerOfMass = this->GetCenterOfMass(it_region, m_InvalidRegions[smallRegions[i]].second, false);
1075  comOfRegions.push_back(centerOfMass);
1076  }
1077 
1078  //Calculate the distance between the center of mass of our small region and all other regions
1079  itk::ImageRegionIterator<IntegerImageType> it_regionSmallest(m_Image, m_InvalidRegions[smallestRegions[thisSmallRegion]].first);
1080  std::vector<double> centerOfMassCurrentRegion = this->GetCenterOfMass(it_regionSmallest, m_InvalidRegions[smallestRegions[thisSmallRegion]].second, false);
1081 
1082  distance.clear();
1083 
1084  for (int i = 0; i < smallRegions.size(); i++)
1085  {
1086  itk::ImageRegionIterator<IntegerImageType> it_region(m_Image, m_InvalidRegions[smallRegions[i]].first);
1087  centerOfMassPossibleRegion = this->GetCenterOfMass(it_region, m_InvalidRegions[smallRegions[i]].second, false);
1088  distance.push_back(this->GetDistanceVector(centerOfMassCurrentRegion, centerOfMassPossibleRegion));
1089  }
1090  smallDistances.clear();
1091  std::vector<double> distanceWithoutZero = distance;
1092  for (int i = 0; i < distanceWithoutZero.size(); i++)
1093  {
1094  if (distanceWithoutZero[i] == 0)
1095  {
1096  distanceWithoutZero.erase(distanceWithoutZero.begin() + i); //Our smallestRegion is calculated too, erase the zero entry
1097  break;
1098  }
1099  }
1100 
1101  if (distanceWithoutZero.size() > 0)
1102  {
1103  //If there are small Regions with equal distance btw the COM's choose a random one
1104  for (int i = 0; i < distance.size(); i++)
1105  {
1106  if (distance[i] == this->SmallestValue(distanceWithoutZero))
1107  {
1108  smallDistances.push_back(smallRegions[i]);
1109  }
1110  }
1111 
1112  thisIndex = smallDistances[rand() % smallDistances.size()];
1113 
1114  //Colour the small parcel like the found one
1115  for (it_regionSmallest.GoToBegin(); !it_regionSmallest.IsAtEnd(); ++it_regionSmallest)
1116  {
1117  if (it_regionSmallest.Value() == m_InvalidRegions[smallestRegions[thisSmallRegion]].second)
1118  {
1119  it_regionSmallest.Value() = m_InvalidRegions[thisIndex].second;
1120  }
1121  }
1122 
1123  //Expand the region of the new parcel if necessary!
1124  chosenRegion = m_InvalidRegions[thisIndex];
1125  chosenRegion.first = this->ExtendedRegionNotNeighbor(chosenRegion.first , m_InvalidRegions[smallestRegions[thisSmallRegion]].first);
1126  m_InvalidRegions[thisIndex] = chosenRegion;
1127 
1128  //Erase the smallest Region from m_InvalidRegions
1129  m_InvalidRegions.erase(m_InvalidRegions.begin() + smallestRegions[thisSmallRegion]);
1130  currentNumberOfParcels--;
1131 
1132  if (currentNumberOfParcels <= m_DesiredNumberOfParcels)
1133  {
1134  tooManyParcels = false;
1135  }
1136 
1137  //m_SizeOfRegions changed, get the new one
1138  this->GetSizeOfRegions();
1139  }
1140 
1141  else//No merging possible, erase this small region from m_InvalidRegions and from m_SizeOfRegions
1142  {
1143  m_InvalidRegions.erase(m_InvalidRegions.begin() + smallestRegions[thisSmallRegion]);
1144  currentNumberOfParcels--;
1145  m_SizeOfFinishedRegions.push_back(m_SizeOfRegions[smallestRegions[thisSmallRegion]]);
1146  m_SizeOfRegions.erase(m_SizeOfRegions.begin() + smallestRegions[thisSmallRegion]);
1147 
1148  if (currentNumberOfParcels <= m_DesiredNumberOfParcels)
1149  {
1150  tooManyParcels = false;
1151  }
1152  }
1153  }
1154  }
1155  else //there are no parcels left to merge with
1156  {
1157  mergingStillPossible = false;
1158  }
1159  }
1160  return (m_SizeOfFinishedRegions.size() + m_SizeOfRegions.size());
1161 }
1162 
1163 template < typename TPixel, unsigned int VImageDimension >
1165 {
1166  m_SizeOfRegions.clear();
1168  voxelCounter.SetImage(m_Image);
1169  for (int i = 0; i < m_InvalidRegions.size(); i++)
1170  {
1171  voxelCounter.SetRegion(m_InvalidRegions[i].first);
1172  m_SizeOfRegions.push_back(voxelCounter.VoxelWithValue(m_InvalidRegions[i].second));
1173  }
1174 }
1175 
1176 template < typename TPixel, unsigned int VImageDimension >
1178 {
1180  voxelCounter.SetImage(m_Image);
1181  std::stringstream message;
1182  for (int i = 0; i < m_InvalidRegions.size(); i++)
1183  {
1184  voxelCounter.SetRegion(m_InvalidRegions[i].first);
1185  m_SizeOfRegions.push_back(voxelCounter.VoxelWithValue(m_InvalidRegions[i].second));
1186  message << voxelCounter.VoxelWithValue(m_InvalidRegions[i].second) << " , ";
1187  }
1188  MITK_DEBUG << message.str();
1189 }
void SetImage(itk::Image< TPixel, VImageDimension > *)
void SetRegion(std::pair< RegionType, int >)
void SetImage(itk::Image< TPixel, VImageDimension > *)
int MergeParcels()
Merge parcels according to a cost function Looks for the parcel with the smallest number of voxels...
itk::Image< unsigned char, 3 > ImageType
bool IsUnique(int number, std::vector< int > vec)
Checks if a number is an element of the vector already.
#define MITK_DEBUG
Definition: mitkLogMacros.h:26
double GetDistance(std::vector< double > centerOfMass, typename ImageType::IndexType indexNewVoxel)
Calculates the distance between two voxels, the position of the first one is given by an index and th...
int CalculateCost()
Checks if the functions RegionIsFullEnough, SidesAreNotTooLong and COMLiesWithinParcel all return val...
int MaximalValue()
Calculates the maximal possible value of the cost function.
static Vector3D offset
void SetImage(itk::Image< TPixel, VImageDimension > *)
void ShowSizeOfRegions()
Calculates and shows the size (number of voxels) of all regions on the console.
void SetRegion(typename ImageType::RegionType)
unsigned int xrand(void)
void SetBoolsForMerging(bool mergingWithNumberParcels, bool mergingWithSmallestParcel, bool justMergeSmallParcels)
void SetVariablesForMerging(int givenSizeOfSmallestRegion, int desiredNumberOfParcels, int givenSizeOfSmallestRegionBeginning)
std::vector< double > GetCenterOfMass(itk::ImageRegionIterator< ImageType > it_region, int valueOfRegion, bool getSizeOfRegions)
Gives back the center of mass and -if wanted- the size (number of voxels) of a parcel.
int VoxelWithValue(int value)
Counts all voxels with the chosen value in the set region.
static T max(T x, T y)
Definition: svm.cpp:70
void GetSizeOfRegions()
Calculates the size (number of voxels) of all regions.
void SetAppropriateValues()
Changes the values of the nodes, such that no gaps exist and it starts with value 1...
int SmallestValue(std::vector< int > distance)
Gives back the smallest value of an int-vector.
ImageType::RegionType ExtendedRegion(typename ImageType::RegionType chosenRegion, typename ImageType::IndexType indexChosenVoxel)
Extends the region if the chosen voxel lies outside.
static T min(T x, T y)
Definition: svm.cpp:67
void GetRandomSeedVoxels()
Sets randomly chosen seed voxels (1x1x1 regions) on the segmented image This is done by creating a ve...
double GetDistanceVector(std::vector< double > centerOfMass, std::vector< double > indexNewVoxel)
Calculates the distance between two voxels, both positions are given by vectors.
void AllocateIsolatedVoxels()
Add voxels of the segmented part to an appropriate region (no neighbors necessary) Checks which voxel...
ImageType::RegionType ExtendedRegionNotNeighbor(typename ImageType::RegionType chosenRegion, typename ImageType::RegionType smallestRegion)
Extends the region of a parcel such that the second region lies within.
void FillOverEdgeOrVertex()
Add voxels of the segmented part to an appropriate region (26-connected neighborhood) Checks which vo...
void ParcelGrowthOverFaces()
Add appropriate voxels of the segmented part to a region (just 6-connected neighborhood) A voxel is a...