Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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.