Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkConnectednessFilter.cxx
Go to the documentation of this file.
1 #pragma warning(disable : 4996)
2 
3 #ifndef __ConnectednessFilter_hxx
4 #define __ConnectednessFilter_hxx
5 
7 #include "itkImageRegionConstIterator.h"
8 #include "itkImageRegionIterator.h"
9 #include "itkObjectFactory.h"
10 #include <itkImageRegionConstIteratorWithIndex.h>
11 // Vector Fields
12 #include "itkVectorImage.h"
13 #include <vnl/vnl_vector.h>
14 //
15 #include "itkConstantBoundaryCondition.h"
16 #include "itkNeighborhoodIterator.h"
17 
18 namespace itk
19 {
20  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
22  {
23  this->SetNthInput(0, const_cast<TFeatureImage *>(image));
24  }
25 
26  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
28  {
29  this->SetNthInput(1, const_cast<TSeedImage *>(image));
30  }
31 
32  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
34  {
35  this->SetNthInput(2, const_cast<TSeedImage *>(image));
36  }
37 
38  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
40  const TFeatureImage *conf)
41  {
42  // DO NOT set it as nth input, since it is dti derived and therefore does
43  // not necessarily occupy the same space as the other images
44  m_ConfidenceImage = const_cast<TFeatureImage *>(conf);
45  }
46 
47  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
49  {
50  this->SetNthInput(3, const_cast<TFeatureImage *>(image));
51  }
52 
53  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
55  {
56  return m_DistanceImage;
57  }
58 
59  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
60  typename TFeatureImage::Pointer
62  {
63  return m_EuclideanDistance;
64  }
65 
66  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
68  {
69  // DO NOT set it as nth input, since it is dti derived and therefore does
70  // not necessary occupy the same space as the other images
71  m_TensorImage = const_cast<TTensorImage *>(vecs);
72  }
73 
74  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
76  {
77  const char untouchedValue = -10;
78  const char inListValue = -30;
79  // Get Inputs
80  m_InputImage = dynamic_cast<TFeatureImage *>(this->ProcessObject::GetInput(0));
81  typename TSeedImage::Pointer seed = dynamic_cast<TSeedImage *>(this->ProcessObject::GetInput(1));
82  typename TSeedImage::Pointer mask = dynamic_cast<TSeedImage *>(this->ProcessObject::GetInput(2));
83  typename TFeatureImage::Pointer propagationImage = dynamic_cast<TFeatureImage *>(this->ProcessObject::GetInput(3));
84 
85  if (propagationImage.IsNull())
86  propagationImage = m_InputImage;
87 
88  // ---
89  // Allocate Output
90  typename TFeatureImage::Pointer output = this->GetOutput();
91  output->SetRegions(m_InputImage->GetLargestPossibleRegion());
92  output->Allocate();
93  output->FillBuffer(untouchedValue);
94 
95  m_DistanceImage = TFeatureImage::New();
96  m_DistanceImage->SetRegions(m_InputImage->GetLargestPossibleRegion());
97  m_DistanceImage->Allocate();
98  m_DistanceImage->FillBuffer(untouchedValue);
99  m_DistanceImage->SetOrigin(m_InputImage->GetOrigin());
100  m_DistanceImage->SetDirection(m_InputImage->GetDirection());
101  m_DistanceImage->SetSpacing(m_InputImage->GetSpacing());
102 
103  m_EuclideanDistance = TFeatureImage::New();
104  m_EuclideanDistance->SetRegions(m_InputImage->GetLargestPossibleRegion());
105  m_EuclideanDistance->Allocate();
106  m_EuclideanDistance->FillBuffer(0.0);
107  m_EuclideanDistance->SetOrigin(m_InputImage->GetOrigin());
108  m_EuclideanDistance->SetDirection(m_InputImage->GetDirection());
109  m_EuclideanDistance->SetSpacing(m_InputImage->GetSpacing());
110 
111  // Helper Object - Indicates which voxels are currently in the FIFO,
112  // this prevents excessive memory and comptutational overhead
113  typename TFeatureImage::Pointer activeSeeds = TFeatureImage::New();
114  activeSeeds->SetRegions(m_InputImage->GetLargestPossibleRegion());
115  activeSeeds->Allocate();
116  activeSeeds->FillBuffer(untouchedValue);
117 
118  // Create Iterators
119  SeedIteratorType seedIt(seed, seed->GetLargestPossibleRegion());
120  FeatureIteratorType outIt(output, output->GetLargestPossibleRegion());
121  FeatureIteratorType inputIt(m_InputImage, m_InputImage->GetLargestPossibleRegion());
122 
123  // Set up Neighborhood Iterator to use a 3x3x3 neighborhood
124  typename itk::NeighborhoodIterator<TSeedImage>::RadiusType radius;
125  radius.Fill(1);
126  radius[2] = 1;
127 
128  typename itk::NeighborhoodIterator<TSeedImage> neighbIt(radius, seed, seed->GetRequestedRegion());
129  typename itk::NeighborhoodIterator<TFeatureImage> medianIt(
130  radius, propagationImage, propagationImage->GetRequestedRegion());
131 
132  // Copy Input Image values from SeedRegion values into Output Image
133  // Collect border voxel indices to initialize seed list
134  typedef std::list<IndexType> ContainerType;
135  ContainerType fifo;
136  while (!seedIt.IsAtEnd())
137  {
138  if (seedIt.Get() != 0)
139  {
140  outIt.Set(inputIt.Get());
141  if (activeSeeds->GetPixel(seedIt.GetIndex()) == untouchedValue)
142  {
143  neighbIt.SetLocation(seedIt.GetIndex());
144  for (unsigned int i = 0; i < neighbIt.Size(); ++i)
145  {
146  bool isInside = true;
147  neighbIt.GetPixel(i, isInside);
148  if (!isInside)
149  continue;
150  if (neighbIt.GetPixel(i) == 0) // border voxel
151  {
152  // add index of border voxel into FIFO
153  fifo.push_back(seedIt.GetIndex());
154  activeSeeds->SetPixel(seedIt.GetIndex(), inListValue);
155  // write propagation value into output image
156 
157  if (m_ApplyRankFilter)
158  {
159  // Calculate median in 3x3x3 neighborhood and use this value to propagate
160  medianIt.SetLocation(seedIt.GetIndex());
161  std::vector<FeaturePixelType> samples;
162  for (unsigned int j = 0; j < medianIt.Size(); ++j)
163  {
164  if (seed->GetPixel(medianIt.GetIndex(j)) != 0) // inside seed region
165  samples.push_back(medianIt.GetPixel(j));
166  }
167  std::sort(samples.begin(), samples.end());
168 
169  output->SetPixel(seedIt.GetIndex(), samples.at(samples.size() / 2));
170  }
171  else
172  output->SetPixel(seedIt.GetIndex(), propagationImage->GetPixel(seedIt.GetIndex()));
173 
174  m_DistanceImage->SetPixel(seedIt.GetIndex(), 0);
175  break;
176  }
177  }
178  }
179  }
180 
181  ++inputIt;
182  ++outIt;
183  ++seedIt;
184  }
185  //
186  // Now iterate over items in fifo and check all possible paths starting from these indices
187  // to their neighbors and mark those with lowest cost
188  while (!fifo.empty())
189  {
190  IndexType index = fifo.front();
191  fifo.pop_front();
192  // mark voxel as no longer present in fifo
193  activeSeeds->SetPixel(index, untouchedValue);
194 
195  if (mask->GetPixel(index) == 0)
196  continue;
197 
198  neighbIt.SetLocation(index);
199 
200  FeaturePixelType currDistance = m_DistanceImage->GetPixel(index);
201  FeaturePixelType currEuclidDistance = m_EuclideanDistance->GetPixel(index);
202  FeaturePixelType propagateValue = output->GetPixel(index);
203 
204  // compute cost to get to each neighborhood voxel
205  // if we find a lower way to a neighbor voxel, that voxel gets appended to the fifo
206  for (unsigned int i = 0; i < neighbIt.Size(); ++i)
207  {
208  bool isInside = true;
209  neighbIt.GetPixel(i, isInside);
210  if (!isInside || mask->GetPixel(neighbIt.GetIndex(i)) == 0) // || seed->GetPixel(neighbIt.GetIndex(i)) != 0)
211  continue;
212 
213  if (i == neighbIt.Size() / 2) // skip center voxel
214  continue;
215 
216  FeaturePixelType cost = GetDistanceValue(index, neighbIt.GetIndex(i));
217  if (cost == -1)
218  continue;
219  if (currDistance + cost < m_DistanceImage->GetPixel(neighbIt.GetIndex(i)) ||
220  m_DistanceImage->GetPixel(neighbIt.GetIndex(i)) == untouchedValue)
221  {
222  m_DistanceImage->SetPixel(neighbIt.GetIndex(i), currDistance + cost);
223 
224  // Compute Euclidean distance between indices
225  FeaturePixelType eDist =
226  sqrt(std::pow(index[0] - neighbIt.GetIndex(i)[0], 2) * m_InputImage->GetSpacing()[0] +
227  std::pow(index[1] - neighbIt.GetIndex(i)[1], 2) * m_InputImage->GetSpacing()[1] +
228  std::pow(index[2] - neighbIt.GetIndex(i)[2], 2) * m_InputImage->GetSpacing()[2]);
229 
230  m_EuclideanDistance->SetPixel(neighbIt.GetIndex(i), currEuclidDistance + eDist);
231 
232  output->SetPixel(neighbIt.GetIndex(i), propagateValue);
233  if (activeSeeds->GetPixel(neighbIt.GetIndex(i)) == untouchedValue)
234  {
235  fifo.push_back(neighbIt.GetIndex(i));
236  activeSeeds->SetPixel(neighbIt.GetIndex(i), inListValue);
237  }
238  }
239  }
240  }
241  }
242 
243  template <typename TFeatureImage, typename TSeedImage, typename TTensorImage>
245  IndexType idxEnd)
246  {
247  FeaturePixelType cost = 0;
248 
249  switch (m_Mode)
250  {
251  case (FeatureSimilarity):
252  {
253  cost = (-1 + exp(std::fabs(m_InputImage->GetPixel(idxStart) - m_InputImage->GetPixel(idxEnd))));
254  break;
255  }
256  case (VectorAgreement):
257  {
258  // Rational
259  // Evaluate Confidence of Tensors (startIdx and endIdx), the higher the lower the cost
260  // Evaluate Path direction correspondence with tensor at start and end index
261  vnl_vector_fixed<double, 3> projection;
262  vnl_vector_fixed<double, 3> projectionEnd;
263  projection.fill(0);
264  projectionEnd.fill(0);
265 
266  itk::Point<double, 3> worldPosStart;
267  itk::Point<double, 3> worldPosEnd;
268  IndexTensorType tensorIdx;
269  IndexType featureIdxStart;
270  IndexTensorType tensorIdxEnd;
271  IndexType featureIdxEnd;
272  m_InputImage->TransformIndexToPhysicalPoint(idxStart, worldPosStart);
273  m_InputImage->TransformIndexToPhysicalPoint(idxEnd, worldPosEnd);
274  m_TensorImage->TransformPhysicalPointToIndex(worldPosStart, tensorIdx);
275  m_TensorImage->TransformPhysicalPointToIndex(worldPosEnd, tensorIdxEnd);
276 
277  m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosStart, featureIdxStart);
278  m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosEnd, featureIdxEnd);
279 
280  if (!m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdxEnd) ||
281  !m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdx))
282  return -1;
283 
284  if (m_ConfidenceImage->GetPixel(featureIdxStart) < .1)
285  return -1;
286 
287  vnl_vector_fixed<double, 3> direction;
288  direction[0] = (-worldPosStart[0] + worldPosEnd[0]);
289  direction[1] = (-worldPosStart[1] + worldPosEnd[1]);
290  direction[2] = (-worldPosStart[2] + worldPosEnd[2]);
291 
292  direction = direction.normalize();
293 
294  TensorPixelType tensorEnd = m_TensorImage->GetPixel(tensorIdxEnd);
295  TensorPixelType tensorStart = m_TensorImage->GetPixel(tensorIdx);
296 
297  typename TensorPixelType::EigenValuesArrayType eigenvalues;
298 
299  tensorEnd.ComputeEigenValues(eigenvalues);
300  double maxEVEnd = std::fabs(eigenvalues[2]);
301 
302  tensorStart.ComputeEigenValues(eigenvalues);
303  double maxEVStart = std::fabs(eigenvalues[2]);
304 
305  if (maxEVEnd == 0.0)
306  maxEVEnd = 1; // no change then ..
307  // Normalize by largest EV
308  tensorEnd[0] /= maxEVEnd;
309  tensorEnd[1] /= maxEVEnd;
310  tensorEnd[2] /= maxEVEnd;
311  tensorEnd[3] /= maxEVEnd;
312  tensorEnd[4] /= maxEVEnd;
313  tensorEnd[5] /= maxEVEnd;
314 
315  if (maxEVStart == 0.0)
316  maxEVStart = 1; // no change then ..
317  // Normalize by largest EV
318  tensorStart[0] /= maxEVStart;
319  tensorStart[1] /= maxEVStart;
320  tensorStart[2] /= maxEVStart;
321  tensorStart[3] /= maxEVStart;
322  tensorStart[4] /= maxEVStart;
323  tensorStart[5] /= maxEVStart;
324 
325  /* Matrix Style
326  * | 0 1 2 |
327  * | X 3 4 |
328  * | X X 5 |
329  */
330  // Multiply vector with tensor, and then take the magnitude of it.
331  projection[0] = direction[0] * tensorStart[0] + direction[1] * tensorStart[1] + direction[2] * tensorStart[2];
332  projection[1] = direction[0] * tensorStart[1] + direction[1] * tensorStart[3] + direction[2] * tensorStart[4];
333  projection[2] = direction[0] * tensorStart[2] + direction[1] * tensorStart[4] + direction[2] * tensorStart[5];
334 
335  projectionEnd[0] = direction[0] * tensorEnd[0] + direction[1] * tensorEnd[1] + direction[2] * tensorEnd[2];
336  projectionEnd[1] = direction[0] * tensorEnd[1] + direction[1] * tensorEnd[3] + direction[2] * tensorEnd[4];
337  projectionEnd[2] = direction[0] * tensorEnd[2] + direction[1] * tensorEnd[4] + direction[2] * tensorEnd[5];
338 
339  if (direction.magnitude() == 0.0 || projectionEnd.magnitude() == 0.0)
340  return -1;
341 
342  // Confidences higher than 1 are considered to be noise (due to the behavior of FA values derived from tensors
343  // with negative eigenvalues)
344  FeaturePixelType confStart = m_ConfidenceImage->GetPixel(featureIdxStart);
345  if (confStart > 1)
346  confStart = .1;
347  FeaturePixelType confEnd = m_ConfidenceImage->GetPixel(featureIdxEnd);
348  if (confEnd > 1)
349  confEnd = .1;
350 
351  // costs
352  cost = (1.1 - confStart) * (1.0 / sqrt(projection.magnitude()));
353  cost *= (1.1 - confEnd) * (1.0 / sqrt(projectionEnd.magnitude()));
354 
355  break;
356  }
357  case (FeatureVectorAgreement):
358  {
359  // Rational
360  // Evaluate Confidence of Tensors (startIdx and endIdx), the higher the lower the cost
361  // Evaluate Path direction correspondence with tensor at start and end index
362  vnl_vector_fixed<double, 3> projection;
363  vnl_vector_fixed<double, 3> projectionEnd;
364  projection.fill(0);
365  projectionEnd.fill(0);
366 
367  itk::Point<double, 3> worldPosStart;
368  itk::Point<double, 3> worldPosEnd;
369  IndexTensorType tensorIdx;
370  IndexType featureIdxStart;
371  IndexTensorType tensorIdxEnd;
372  IndexType featureIdxEnd;
373  m_InputImage->TransformIndexToPhysicalPoint(idxStart, worldPosStart);
374  m_InputImage->TransformIndexToPhysicalPoint(idxEnd, worldPosEnd);
375  m_TensorImage->TransformPhysicalPointToIndex(worldPosStart, tensorIdx);
376  m_TensorImage->TransformPhysicalPointToIndex(worldPosEnd, tensorIdxEnd);
377 
378  m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosStart, featureIdxStart);
379  m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosEnd, featureIdxEnd);
380 
381  if (!m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdxEnd) ||
382  !m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdx))
383  return -1;
384 
385  if (m_ConfidenceImage->GetPixel(featureIdxStart) < .1)
386  return -1;
387 
388  vnl_vector_fixed<double, 3> direction;
389  direction[0] = (-worldPosStart[0] + worldPosEnd[0]);
390  direction[1] = (-worldPosStart[1] + worldPosEnd[1]);
391  direction[2] = (-worldPosStart[2] + worldPosEnd[2]);
392 
393  direction = direction.normalize();
394 
395  TensorPixelType tensorEnd = m_TensorImage->GetPixel(tensorIdxEnd);
396  TensorPixelType tensorStart = m_TensorImage->GetPixel(tensorIdx);
397 
398  typename TensorPixelType::EigenValuesArrayType eigenvalues;
399 
400  tensorEnd.ComputeEigenValues(eigenvalues);
401  double maxEVEnd = std::fabs(eigenvalues[2]);
402 
403  tensorStart.ComputeEigenValues(eigenvalues);
404  double maxEVStart = std::fabs(eigenvalues[2]);
405 
406  if (maxEVEnd == 0.0)
407  maxEVEnd = 1; // no change then ..
408  // Normalize by largest EV
409  tensorEnd[0] /= maxEVEnd;
410  tensorEnd[1] /= maxEVEnd;
411  tensorEnd[2] /= maxEVEnd;
412  tensorEnd[3] /= maxEVEnd;
413  tensorEnd[4] /= maxEVEnd;
414  tensorEnd[5] /= maxEVEnd;
415 
416  if (maxEVStart == 0.0)
417  maxEVStart = 1; // no change then ..
418  // Normalize by largest EV
419  tensorStart[0] /= maxEVStart;
420  tensorStart[1] /= maxEVStart;
421  tensorStart[2] /= maxEVStart;
422  tensorStart[3] /= maxEVStart;
423  tensorStart[4] /= maxEVStart;
424  tensorStart[5] /= maxEVStart;
425 
426  /* Matrix Style
427  * | 0 1 2 |
428  * | X 3 4 |
429  * | X X 5 |
430  */
431  // Multiply vector with tensor, and then take the magnitude of it.
432  projection[0] = direction[0] * tensorStart[0] + direction[1] * tensorStart[1] + direction[2] * tensorStart[2];
433  projection[1] = direction[0] * tensorStart[1] + direction[1] * tensorStart[3] + direction[2] * tensorStart[4];
434  projection[2] = direction[0] * tensorStart[2] + direction[1] * tensorStart[4] + direction[2] * tensorStart[5];
435 
436  projectionEnd[0] = direction[0] * tensorEnd[0] + direction[1] * tensorEnd[1] + direction[2] * tensorEnd[2];
437  projectionEnd[1] = direction[0] * tensorEnd[1] + direction[1] * tensorEnd[3] + direction[2] * tensorEnd[4];
438  projectionEnd[2] = direction[0] * tensorEnd[2] + direction[1] * tensorEnd[4] + direction[2] * tensorEnd[5];
439 
440  if (direction.magnitude() == 0.0 || projectionEnd.magnitude() == 0.0)
441  return -1;
442 
443  // Confidences higher than 1 are considered to be noise (due to the behavior of FA values derived from tensors
444  // with negative eigenvalues)
445  FeaturePixelType confStart = m_ConfidenceImage->GetPixel(featureIdxStart);
446  if (confStart > 1)
447  confStart = .1;
448  FeaturePixelType confEnd = m_ConfidenceImage->GetPixel(featureIdxEnd);
449  if (confEnd > 1)
450  confEnd = .1;
451 
452  // costs
453  cost = (1.1 - confStart) * (1.0 / sqrt(projection.magnitude()));
454  cost *= (1.1 - confEnd) * (1.0 / sqrt(projectionEnd.magnitude()));
455  cost *= (-1 + exp(std::fabs(m_InputImage->GetPixel(idxStart) - m_InputImage->GetPixel(idxEnd))));
456 
457  break;
458  }
459  }
460  return cost;
461  }
462 
463 } // end itk namespace
464 
465 #endif
itk::ImageRegionIterator< TSeedImage > SeedIteratorType
itk::SmartPointer< Self > Pointer
void SetInputVectorField(const TTensorImage *vecs)
SetInputVectorField.
TFeatureImage::IndexType IndexType
TTensorImage::PixelType TensorPixelType
void SetInputSeed(const TSeedImage *mask)
SetInputSeed - seed area from which features are to be propagated.
itk::Image< itk::DiffusionTensor3D< TTensorImagePixelType >, 3 > TTensorImage
void SetInputImage(const TFeatureImage *image)
itk::ImageRegionIterator< TFeatureImage > FeatureIteratorType
void SetInputMask(const TSeedImage *mask)
SetInputMask - Filter only operates in masked area.
void SetPropagationImage(const TFeatureImage *prop)
SetPropagationImage - OPTIONAL. Set image which values are propagated, it this is not supplied the In...
TTensorImage::IndexType IndexTensorType
void SetInputVectorFieldConfidenceMap(const TFeatureImage *conf)
SetInputVectorFieldConfidenceMap - Map that assigned a weight/confidence to each vector.
TFeatureImage::Pointer GetDistanceImage()
GetDistanceImage - Return the distance image using the specified metric.
TFeatureImage::PixelType FeaturePixelType
virtual double GetDistanceValue(IndexType idxStart, IndexType idxEnd)
TFeatureImage::Pointer GetEuclideanDistanceImage()
GetEuclideanDistanceImage - Return distance image that provides distance of shortest path for each vo...
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.