1 #pragma warning(disable : 4996)
3 #ifndef __ConnectednessFilter_hxx
4 #define __ConnectednessFilter_hxx
7 #include "itkImageRegionConstIterator.h"
8 #include "itkImageRegionIterator.h"
9 #include "itkObjectFactory.h"
10 #include <itkImageRegionConstIteratorWithIndex.h>
12 #include "itkVectorImage.h"
13 #include <vnl/vnl_vector.h>
15 #include "itkConstantBoundaryCondition.h"
16 #include "itkNeighborhoodIterator.h"
20 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
23 this->SetNthInput(0, const_cast<TFeatureImage *>(image));
26 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
29 this->SetNthInput(1, const_cast<TSeedImage *>(image));
32 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
35 this->SetNthInput(2, const_cast<TSeedImage *>(image));
38 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
40 const TFeatureImage *conf)
44 m_ConfidenceImage =
const_cast<TFeatureImage *
>(conf);
47 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
50 this->SetNthInput(3, const_cast<TFeatureImage *>(image));
53 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
56 return m_DistanceImage;
59 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
63 return m_EuclideanDistance;
66 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
74 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
77 const char untouchedValue = -10;
78 const char inListValue = -30;
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));
85 if (propagationImage.IsNull())
86 propagationImage = m_InputImage;
91 output->SetRegions(m_InputImage->GetLargestPossibleRegion());
93 output->FillBuffer(untouchedValue);
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());
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());
114 activeSeeds->SetRegions(m_InputImage->GetLargestPossibleRegion());
115 activeSeeds->Allocate();
116 activeSeeds->FillBuffer(untouchedValue);
124 typename itk::NeighborhoodIterator<TSeedImage>::RadiusType radius;
128 typename itk::NeighborhoodIterator<TSeedImage> neighbIt(radius, seed, seed->GetRequestedRegion());
129 typename itk::NeighborhoodIterator<TFeatureImage> medianIt(
130 radius, propagationImage, propagationImage->GetRequestedRegion());
134 typedef std::list<IndexType> ContainerType;
136 while (!seedIt.IsAtEnd())
138 if (seedIt.Get() != 0)
140 outIt.Set(inputIt.Get());
141 if (activeSeeds->GetPixel(seedIt.GetIndex()) == untouchedValue)
143 neighbIt.SetLocation(seedIt.GetIndex());
144 for (
unsigned int i = 0; i < neighbIt.Size(); ++i)
146 bool isInside =
true;
147 neighbIt.GetPixel(i, isInside);
150 if (neighbIt.GetPixel(i) == 0)
153 fifo.push_back(seedIt.GetIndex());
154 activeSeeds->SetPixel(seedIt.GetIndex(), inListValue);
157 if (m_ApplyRankFilter)
160 medianIt.SetLocation(seedIt.GetIndex());
161 std::vector<FeaturePixelType> samples;
162 for (
unsigned int j = 0; j < medianIt.Size(); ++j)
164 if (seed->GetPixel(medianIt.GetIndex(j)) != 0)
165 samples.push_back(medianIt.GetPixel(j));
167 std::sort(samples.begin(), samples.end());
169 output->SetPixel(seedIt.GetIndex(), samples.at(samples.size() / 2));
172 output->SetPixel(seedIt.GetIndex(), propagationImage->GetPixel(seedIt.GetIndex()));
174 m_DistanceImage->SetPixel(seedIt.GetIndex(), 0);
188 while (!fifo.empty())
193 activeSeeds->SetPixel(index, untouchedValue);
195 if (mask->GetPixel(index) == 0)
198 neighbIt.SetLocation(index);
206 for (
unsigned int i = 0; i < neighbIt.Size(); ++i)
208 bool isInside =
true;
209 neighbIt.GetPixel(i, isInside);
210 if (!isInside || mask->GetPixel(neighbIt.GetIndex(i)) == 0)
213 if (i == neighbIt.Size() / 2)
219 if (currDistance + cost < m_DistanceImage->GetPixel(neighbIt.GetIndex(i)) ||
220 m_DistanceImage->GetPixel(neighbIt.GetIndex(i)) == untouchedValue)
222 m_DistanceImage->SetPixel(neighbIt.GetIndex(i), currDistance + cost);
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]);
230 m_EuclideanDistance->SetPixel(neighbIt.GetIndex(i), currEuclidDistance + eDist);
232 output->SetPixel(neighbIt.GetIndex(i), propagateValue);
233 if (activeSeeds->GetPixel(neighbIt.GetIndex(i)) == untouchedValue)
235 fifo.push_back(neighbIt.GetIndex(i));
236 activeSeeds->SetPixel(neighbIt.GetIndex(i), inListValue);
243 template <
typename TFeatureImage,
typename TSeedImage,
typename TTensorImage>
251 case (FeatureSimilarity):
253 cost = (-1 + exp(std::fabs(m_InputImage->GetPixel(idxStart) - m_InputImage->GetPixel(idxEnd))));
256 case (VectorAgreement):
261 vnl_vector_fixed<double, 3> projection;
262 vnl_vector_fixed<double, 3> projectionEnd;
264 projectionEnd.fill(0);
266 itk::Point<double, 3> worldPosStart;
267 itk::Point<double, 3> worldPosEnd;
272 m_InputImage->TransformIndexToPhysicalPoint(idxStart, worldPosStart);
273 m_InputImage->TransformIndexToPhysicalPoint(idxEnd, worldPosEnd);
274 m_TensorImage->TransformPhysicalPointToIndex(worldPosStart, tensorIdx);
275 m_TensorImage->TransformPhysicalPointToIndex(worldPosEnd, tensorIdxEnd);
277 m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosStart, featureIdxStart);
278 m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosEnd, featureIdxEnd);
280 if (!m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdxEnd) ||
281 !m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdx))
284 if (m_ConfidenceImage->GetPixel(featureIdxStart) < .1)
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]);
292 direction = direction.normalize();
297 typename TensorPixelType::EigenValuesArrayType eigenvalues;
299 tensorEnd.ComputeEigenValues(eigenvalues);
300 double maxEVEnd = std::fabs(eigenvalues[2]);
302 tensorStart.ComputeEigenValues(eigenvalues);
303 double maxEVStart = std::fabs(eigenvalues[2]);
308 tensorEnd[0] /= maxEVEnd;
309 tensorEnd[1] /= maxEVEnd;
310 tensorEnd[2] /= maxEVEnd;
311 tensorEnd[3] /= maxEVEnd;
312 tensorEnd[4] /= maxEVEnd;
313 tensorEnd[5] /= maxEVEnd;
315 if (maxEVStart == 0.0)
318 tensorStart[0] /= maxEVStart;
319 tensorStart[1] /= maxEVStart;
320 tensorStart[2] /= maxEVStart;
321 tensorStart[3] /= maxEVStart;
322 tensorStart[4] /= maxEVStart;
323 tensorStart[5] /= maxEVStart;
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];
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];
339 if (direction.magnitude() == 0.0 || projectionEnd.magnitude() == 0.0)
352 cost = (1.1 - confStart) * (1.0 / sqrt(projection.magnitude()));
353 cost *= (1.1 - confEnd) * (1.0 / sqrt(projectionEnd.magnitude()));
357 case (FeatureVectorAgreement):
362 vnl_vector_fixed<double, 3> projection;
363 vnl_vector_fixed<double, 3> projectionEnd;
365 projectionEnd.fill(0);
367 itk::Point<double, 3> worldPosStart;
368 itk::Point<double, 3> worldPosEnd;
373 m_InputImage->TransformIndexToPhysicalPoint(idxStart, worldPosStart);
374 m_InputImage->TransformIndexToPhysicalPoint(idxEnd, worldPosEnd);
375 m_TensorImage->TransformPhysicalPointToIndex(worldPosStart, tensorIdx);
376 m_TensorImage->TransformPhysicalPointToIndex(worldPosEnd, tensorIdxEnd);
378 m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosStart, featureIdxStart);
379 m_ConfidenceImage->TransformPhysicalPointToIndex(worldPosEnd, featureIdxEnd);
381 if (!m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdxEnd) ||
382 !m_TensorImage->GetLargestPossibleRegion().IsInside(tensorIdx))
385 if (m_ConfidenceImage->GetPixel(featureIdxStart) < .1)
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]);
393 direction = direction.normalize();
398 typename TensorPixelType::EigenValuesArrayType eigenvalues;
400 tensorEnd.ComputeEigenValues(eigenvalues);
401 double maxEVEnd = std::fabs(eigenvalues[2]);
403 tensorStart.ComputeEigenValues(eigenvalues);
404 double maxEVStart = std::fabs(eigenvalues[2]);
409 tensorEnd[0] /= maxEVEnd;
410 tensorEnd[1] /= maxEVEnd;
411 tensorEnd[2] /= maxEVEnd;
412 tensorEnd[3] /= maxEVEnd;
413 tensorEnd[4] /= maxEVEnd;
414 tensorEnd[5] /= maxEVEnd;
416 if (maxEVStart == 0.0)
419 tensorStart[0] /= maxEVStart;
420 tensorStart[1] /= maxEVStart;
421 tensorStart[2] /= maxEVStart;
422 tensorStart[3] /= maxEVStart;
423 tensorStart[4] /= maxEVStart;
424 tensorStart[5] /= maxEVStart;
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];
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];
440 if (direction.magnitude() == 0.0 || projectionEnd.magnitude() == 0.0)
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))));
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
virtual void GenerateData()
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.