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