diff --git a/Core/Code/Resources/Interactions/Legacy/StateMachine.xml b/Core/Code/Resources/Interactions/Legacy/StateMachine.xml
index 7ce0de4ab3..0cd6f2be12 100644
--- a/Core/Code/Resources/Interactions/Legacy/StateMachine.xml
+++ b/Core/Code/Resources/Interactions/Legacy/StateMachine.xml
@@ -1,4061 +1,4101 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Modules/CMakeLists.txt b/Modules/CMakeLists.txt
index 1220512bd6..dc1781f1a1 100644
--- a/Modules/CMakeLists.txt
+++ b/Modules/CMakeLists.txt
@@ -1,60 +1,60 @@
set(LIBPOSTFIX "Ext")
# Modules must be listed according to their dependencies
set(module_dirs
SceneSerializationBase
PlanarFigure
ImageExtraction
ImageStatistics
LegacyAdaptors
IpPicSupport
MitkExt
SceneSerialization
+ GraphAlgorithms
Segmentation
Qmitk
QmitkExt
Properties
- GraphAlgorithms
DiffusionImaging
GPGPU
IGT
CameraCalibration
IGTUI
RigidRegistration
RigidRegistrationUI
DeformableRegistration
DeformableRegistrationUI
OpenCL
OpenCVVideoSupport
Overlays
InputDevices
ToFHardware
ToFProcessing
ToFUI
US
ClippingTools
USUI
DicomUI
Simulation
Python
)
set(MITK_DEFAULT_SUBPROJECTS MITK-Modules)
foreach(module_dir ${module_dirs})
add_subdirectory(${module_dir})
endforeach()
if(MITK_PRIVATE_MODULES)
file(GLOB all_subdirs RELATIVE ${MITK_PRIVATE_MODULES} ${MITK_PRIVATE_MODULES}/*)
foreach(subdir ${all_subdirs})
string(FIND ${subdir} "." _result)
if(_result EQUAL -1)
if(EXISTS ${MITK_PRIVATE_MODULES}/${subdir}/CMakeLists.txt)
message(STATUS "Found private module ${subdir}")
add_subdirectory(${MITK_PRIVATE_MODULES}/${subdir} private_modules/${subdir})
endif()
endif()
endforeach()
endif(MITK_PRIVATE_MODULES)
diff --git a/Modules/GraphAlgorithms/files.cmake b/Modules/GraphAlgorithms/files.cmake
index 0de34bc90f..433112a1c8 100644
--- a/Modules/GraphAlgorithms/files.cmake
+++ b/Modules/GraphAlgorithms/files.cmake
@@ -1,9 +1,10 @@
set(CPP_FILES
itkShortestPathNode.cpp
)
set(H_FILES
itkShortestPathCostFunction.h
itkShortestPathCostFunctionTbss.h
itkShortestPathNode.h
itkShortestPathImageFilter.h
+ itkShortestPathCostFunctionLiveWire.h
)
diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h
new file mode 100644
index 0000000000..033c92951c
--- /dev/null
+++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.h
@@ -0,0 +1,204 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef __itkShortestPathCostFunctionLiveWire_h
+#define __itkShortestPathCostFunctionLiveWire_h
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkShortestPathCostFunction.h" // Superclass of Metrics
+
+#include
+
+
+
+
+namespace itk
+{
+ /** \brief Cost function for LiveWire purposes.
+ Specific features are considered to calculate cummulative
+ costs of a link between two pixels. These are:
+
+ - Gradient Magnitude
+ - Gradient Direction
+ - Laplacian Zero Crossing
+
+ By default the Gradient Magnitude is mapped linear to costs
+ between 0 (good) and 1 (bad). Via SetDynamicCostMap( std::map< int, int > &costMap)
+ a cost map can be set to dynamically map Gradient Magnitude (non
+ linear). Thus lower values can be considered with lower costs
+ than higher values of gradient magnitudes.
+ To compute the costs of the gradient magnitude dynamically
+ a iverted map of the histogram of gradient magnitude image is used.
+
+ */
+ template
+ class ITK_EXPORT ShortestPathCostFunctionLiveWire : public ShortestPathCostFunction
+ {
+ public:
+ /** Standard class typedefs. */
+ typedef ShortestPathCostFunctionLiveWire Self;
+ typedef ShortestPathCostFunction Superclass;
+ typedef SmartPointer Pointer;
+ typedef SmartPointer ConstPointer;
+ typedef itk::ImageRegionConstIterator ConstIteratorType;
+
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(ShortestPathCostFunctionLiveWire, ShortestPathCostFunction);
+
+
+ typedef itk::Image UnsignedCharImageType;
+ typedef itk::Image FloatImageType;
+
+ typedef float ComponentType;
+ typedef itk::CovariantVector< ComponentType, 2 > OutputPixelType;
+ typedef itk::Image< OutputPixelType, 2 > VectorOutputImageType;
+
+ typedef typename TInputImageType::IndexType IndexType;
+ typedef TInputImageType ImageType;
+ typedef itk::ImageRegion<2> RegionType;
+
+
+
+ /** \brief calculates the costs for going from p1 to p2*/
+ virtual double GetCost(IndexType p1, IndexType p2);
+
+ /** \brief returns the minimal costs possible (needed for A*)*/
+ virtual double GetMinCost();
+
+ /** \brief Initialize the metric*/
+ virtual void Initialize ();
+
+
+ /** \brief Set repulsive path*/
+ virtual void AddRepulsivePoint( itk::Index<3> );
+
+ /** \brief Clear repulsive path*/
+ virtual void ClearRepulsivePoints( );
+
+
+
+ ShortestPathCostFunctionLiveWire();
+
+
+ itkSetMacro (RequestedRegion, RegionType);
+ itkGetMacro (RequestedRegion, RegionType);
+
+ // Set/Get function for sigma parameter
+ itkSetMacro (UseApproximateGradient, bool);
+ itkGetMacro (UseApproximateGradient, bool);
+
+ virtual void SetImage(const TInputImageType* _arg)
+ {
+ if (this->m_Image != _arg)
+ {
+ this->m_Image = _arg;
+ this->Modified();
+ this->m_Initialized = false;
+ }
+ }
+
+ void SetDynamicCostMap( std::map< int, int > &costMap)
+ {
+ this->m_CostMap = costMap;
+ this->m_UseCostMap = true;
+ this->m_MaxMapCosts = -1;
+ this->Modified();
+ }
+
+ void SetUseCostMap(bool useCostMap)
+ {
+ this->m_UseCostMap = useCostMap;
+ }
+
+ /**
+ \brief Set the maximum of the dynamic cost map to save computation time.
+ */
+ void SetCostMapMaximum(double max)
+ {
+ this->m_MaxMapCosts = max;
+ }
+
+
+ enum Constants{
+ MAPSCALEFACTOR = 10
+ };
+
+ /** \brief Returns the y value of gaussian with given offset and amplitude
+
+ gaussian approximation
+ f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 )
+
+ \param x
+ \param xOfGaussian - offset
+ \param yOfGaussian - amplitude
+ */
+ static double Gaussian(double x, double xOfGaussian, double yOfGaussian);
+
+ protected:
+
+ virtual ~ShortestPathCostFunctionLiveWire() {};
+
+
+ typename ImageType::Pointer m_GradientMagnImage;
+ typename UnsignedCharImageType::Pointer m_ZeroCrossingsImage;
+ typename FloatImageType::Pointer m_EdgeImage;
+ typename VectorOutputImageType::Pointer m_GradientImage;
+
+ double minCosts;
+
+ bool m_UseRepulsivePoint;
+
+ std::vector< itk::Index<3> > m_RepulsivePoints;
+
+ typename Superclass::PixelType val;
+
+ typename Superclass::PixelType startValue;
+ typename Superclass::PixelType endValue;
+
+ double m_GradientMax;
+
+ RegionType m_RequestedRegion;
+
+ bool m_UseApproximateGradient;
+
+ bool m_Initialized;
+
+ std::map< int, int > m_CostMap;
+
+ bool m_UseCostMap;
+
+ double m_MaxMapCosts;
+
+ private:
+
+ double SigmoidFunction(double I, double max, double min, double alpha, double beta);
+
+
+ };
+
+} // end namespace itk
+
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkShortestPathCostFunctionLiveWire.txx"
+#endif
+
+#endif /* __itkShortestPathCostFunctionLiveWire_h */
diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx
new file mode 100644
index 0000000000..a9ad943a75
--- /dev/null
+++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionLiveWire.txx
@@ -0,0 +1,441 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef __itkShortestPathCostFunctionLiveWire_txx
+#define __itkShortestPathCostFunctionLiveWire_txx
+
+#include "itkShortestPathCostFunctionLiveWire.h"
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+namespace itk
+{
+
+ // Constructor
+ template
+ ShortestPathCostFunctionLiveWire
+ ::ShortestPathCostFunctionLiveWire()
+ {
+ m_UseRepulsivePoint = false;
+ m_GradientMax = 0.0;
+ m_Initialized = false;
+ m_UseCostMap = false;
+ m_MaxMapCosts = -1.0;
+ }
+
+
+
+ template
+ void ShortestPathCostFunctionLiveWire
+ ::AddRepulsivePoint( itk::Index<3> c )
+ {
+ m_RepulsivePoints.push_back(c);
+ m_UseRepulsivePoint = true;
+ }
+
+
+
+ template
+ void ShortestPathCostFunctionLiveWire
+ ::ClearRepulsivePoints()
+ {
+ m_RepulsivePoints.clear();
+ m_UseRepulsivePoint = false;
+ }
+
+
+
+ template
+ double ShortestPathCostFunctionLiveWire
+ ::GetCost(IndexType p1 ,IndexType p2)
+ {
+
+ unsigned long xMAX = this->m_Image->GetLargestPossibleRegion().GetSize()[0];
+ unsigned long yMAX = this->m_Image->GetLargestPossibleRegion().GetSize()[1];
+
+ double gradientX, gradientY;
+ gradientX = gradientY = 0.0;
+
+ double gradientCost;
+
+ double gradientMagnitude;
+
+ /* ++++++++++++++++++++ GradientMagnitude costs ++++++++++++++++++++++++++*/
+
+ gradientMagnitude = this->m_GradientMagnImage->GetPixel(p2);
+ gradientX = m_GradientImage->GetPixel(p2)[0];
+ gradientY = m_GradientImage->GetPixel(p2)[1];
+
+
+
+ if(m_UseCostMap && !m_CostMap.empty())
+ {
+ std::map< int, int >::iterator end = m_CostMap.end();
+ std::map< int, int >::iterator last = --(m_CostMap.end());
+
+ //current position
+ std::map< int, int >::iterator x;
+ //std::map< int, int >::key_type keyOfX = static_cast::key_type>(gradientMagnitude * 1000);
+ int keyOfX = static_cast(gradientMagnitude /* ShortestPathCostFunctionLiveWire::MAPSCALEFACTOR*/);
+ x = m_CostMap.find( keyOfX );
+
+ std::map< int, int >::iterator left2;
+ std::map< int, int >::iterator left1;
+ std::map< int, int >::iterator right1;
+ std::map< int, int >::iterator right2;
+
+ if( x == end )
+ {//x can also be == end if the key is not in the map but between two other keys
+ //search next key within map from x upwards
+ right1 = m_CostMap.lower_bound( keyOfX );
+ }
+ else
+ {
+ right1 = x;
+ }
+
+ if(right1 == end || right1 == last )
+ {
+ right2 = end;
+ }
+ else//( right1 != (end-1) )
+ {
+ std::map< int, int >::iterator temp = right1;
+ right2 = ++right1;//rght1 + 1
+ right1 = temp;
+ }
+
+
+ if( right1 == m_CostMap.begin() )
+ {
+ left1 = end;
+ left2 = end;
+ }
+ else if( right1 == (++(m_CostMap.begin())) )
+ {
+ std::map< int, int >::iterator temp = right1;
+ left1 = --right1;//rght1 - 1
+ right1 = temp;
+ left2 = end;
+ }
+ else
+ {
+ std::map< int, int >::iterator temp = right1;
+ left1 = --right1;//rght1 - 1
+ left2 = --right1;//rght1 - 2
+ right1 = temp;
+ }
+
+ double partRight1, partRight2, partLeft1, partLeft2;
+ partRight1 = partRight2 = partLeft1 = partLeft2 = 0.0;
+
+
+ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 )
+
+ gaussian approximation
+
+ where
+ v(bin) is the value in the map
+ k(bin) is the key
+ */
+ if( left2 != end )
+ {
+ partLeft2 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, left2->first, left2->second);
+ }
+
+ if( left1 != end )
+ {
+ partLeft1 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, left1->first, left1->second);
+ }
+
+ if( right1 != end )
+ {
+ partRight1 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, right1->first, right1->second);
+ }
+
+ if( right2 != end )
+ {
+ partRight2 = ShortestPathCostFunctionLiveWire::Gaussian(keyOfX, right2->first, right2->second);
+ }
+ /*----------------------------------------------------------------------------*/
+
+
+ if( m_MaxMapCosts > 0.0 )
+ {
+ gradientCost = 1.0 - ( (partRight1 + partRight2 + partLeft1 + partLeft2) / m_MaxMapCosts );
+ }
+ else
+ {//use linear mapping
+ gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
+ }
+
+ }
+ else
+ {//use linear mapping
+ //value between 0 (good) and 1 (bad)
+ gradientCost = 1.0 - (gradientMagnitude / m_GradientMax);
+
+ }
+ /* -----------------------------------------------------------------------------*/
+
+
+
+ /* ++++++++++++++++++++ Laplacian zero crossing costs ++++++++++++++++++++++++++*/
+ // f(p) = 0; if I(p)=0
+ // or 1; if I(p)!=0
+ double laplacianCost;
+ typename Superclass::PixelType laplaceImageValue;
+
+
+ laplaceImageValue = m_EdgeImage->GetPixel(p2);
+
+ if(laplaceImageValue < 0 || laplaceImageValue > 0)
+ {
+ laplacianCost = 1.0;
+ }
+ else
+ {
+ laplacianCost = 0.0;
+ }
+
+ /* -----------------------------------------------------------------------------*/
+
+
+
+ /* ++++++++++++++++++++ Gradient direction costs ++++++++++++++++++++++++++*/
+ //vector q-p i.e. p2-p1
+ double vQP[2];
+ vQP[0] = p2[0] - p1[0];
+ vQP[1] = p2[1] - p1[1];
+ //-------
+
+ //vector p-q i.e. p1-p2
+ double vPQ[2];
+ vPQ[0] = p1[0] - p2[0];
+ vPQ[1] = p1[1] - p2[1];
+ //-------
+
+ // gradient vector at p1
+ double nGradientAtP1[2];
+ nGradientAtP1[0] = gradientX;//previously computed for gradient magnitude
+ nGradientAtP1[1] = gradientY;
+
+ //gradient direction unit vector of p1
+ nGradientAtP1[0] /= gradientMagnitude;
+ nGradientAtP1[1] /= gradientMagnitude;
+ //-------
+
+ // gradient vector at p1
+ double nGradientAtP2[2];
+
+
+ nGradientAtP2[0] = m_GradientImage->GetPixel(p2)[0];
+ nGradientAtP2[1] = m_GradientImage->GetPixel(p2)[1];
+
+ nGradientAtP2[0] /= m_GradientMagnImage->GetPixel(p2);
+ nGradientAtP2[1] /= m_GradientMagnImage->GetPixel(p2);
+
+
+ double scalarProduct = (nGradientAtP1[0] * nGradientAtP2[0]) + (nGradientAtP1[1] * nGradientAtP2[1]);
+ if( abs(scalarProduct) >= 1.0)
+ {
+ //this should probably not happen; make sure the input for acos is valid
+ scalarProduct = 0.999999999;
+ }
+
+ double gradientDirectionCost = acos( scalarProduct ) / 3.14159265;
+ /*------------------------------------------------------------------------*/
+
+
+
+
+ /*+++++++++++++++++++++ local component costs +++++++++++++++++++++++++++*/
+ /*weights*/
+ double w1;
+ double w2;
+ double w3;
+ double costs = 0.0;
+
+ if (this->m_UseCostMap){
+ w1 = 0.43;
+ w2= 0.43;
+ w3 = 0.14;
+ }else{
+ w1 = 0.10;
+ w2= 0.85;
+ w3 = 0.05;
+ }
+ costs = w1 * laplacianCost + w2 * gradientCost + w3 * gradientDirectionCost;
+
+
+ //scale by euclidian distance
+ double costScale;
+ if( p1[0] == p2[0] || p1[1] == p2[1])
+ {
+ //horizontal or vertical neighbor
+ costScale = 1.0;
+ }
+ else
+ {
+ //diagonal neighbor
+ costScale = sqrt(2.0);
+ }
+
+ costs *= costScale;
+
+ return costs;
+ }
+
+
+
+ template
+ double ShortestPathCostFunctionLiveWire
+ ::GetMinCost()
+ {
+ return minCosts;
+ }
+
+
+
+ template
+ void ShortestPathCostFunctionLiveWire
+ ::Initialize()
+ {
+ if(!m_Initialized)
+ {
+
+
+ typedef itk::CastImageFilter< TInputImageType, FloatImageType > CastFilterType;
+ typename CastFilterType::Pointer castFilter = CastFilterType::New();
+ castFilter->SetInput(this->m_Image);
+
+
+ // init gradient magnitude image
+ typedef itk::GradientMagnitudeImageFilter< FloatImageType, FloatImageType> GradientMagnitudeFilterType;
+ typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New();
+ gradientFilter->SetInput(castFilter->GetOutput());
+ //gradientFilter->SetNumberOfThreads(4);
+ //gradientFilter->GetOutput()->SetRequestedRegion(m_RequestedRegion);
+
+ gradientFilter->Update();
+ m_GradientMagnImage = gradientFilter->GetOutput();
+
+ typedef itk::StatisticsImageFilter StatisticsImageFilterType;
+ typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New();
+ statisticsImageFilter->SetInput(this->m_GradientMagnImage);
+ statisticsImageFilter->Update();
+
+ m_GradientMax = statisticsImageFilter->GetMaximum();
+
+
+
+ //Filter class is instantiated
+ /*typedef itk::GradientRecursiveGaussianImageFilter GradientFilterType;*/
+
+ typedef itk::GradientImageFilter< FloatImageType > GradientFilterType;
+
+ typename GradientFilterType::Pointer filter = GradientFilterType::New();
+ //sigma is specified in millimeters
+ //filter->SetSigma( 1.5 );
+ filter->SetInput(castFilter->GetOutput());
+ filter->Update();
+
+ m_GradientImage = filter->GetOutput();
+
+
+ // init zero crossings
+ //typedef itk::ZeroCrossingImageFilter< TInputImageType, UnsignedCharImageType > ZeroCrossingImageFilterType;
+ //ZeroCrossingImageFilterType::Pointer zeroCrossingImageFilter = ZeroCrossingImageFilterType::New();
+ //zeroCrossingImageFilter->SetInput(this->m_Image);
+ //zeroCrossingImageFilter->SetBackgroundValue(1);
+ //zeroCrossingImageFilter->SetForegroundValue(0);
+ //zeroCrossingImageFilter->SetNumberOfThreads(4);
+ //zeroCrossingImageFilter->Update();
+
+ //m_EdgeImage = zeroCrossingImageFilter->GetOutput();
+
+
+ //cast image to float to apply canny edge dection filter
+ /*typedef itk::CastImageFilter< TInputImageType, FloatImageType > CastFilterType;
+ CastFilterType::Pointer castFilter = CastFilterType::New();
+ castFilter->SetInput(this->m_Image);*/
+
+ //typedef itk::LaplacianImageFilter filterType;
+ //filterType::Pointer laplacianFilter = filterType::New();
+ //laplacianFilter->SetInput( castFilter->GetOutput() ); // NOTE: input image type must be double or float
+ //laplacianFilter->Update();
+
+ //m_EdgeImage = laplacianFilter->GetOutput();
+
+ //init canny edge detection
+ typedef itk::CannyEdgeDetectionImageFilter CannyEdgeDetectionImageFilterType;
+ typename CannyEdgeDetectionImageFilterType::Pointer cannyEdgeDetectionfilter = CannyEdgeDetectionImageFilterType::New();
+ cannyEdgeDetectionfilter->SetInput(castFilter->GetOutput());
+ cannyEdgeDetectionfilter->SetUpperThreshold(30);
+ cannyEdgeDetectionfilter->SetLowerThreshold(15);
+ cannyEdgeDetectionfilter->SetVariance(4);
+ cannyEdgeDetectionfilter->SetMaximumError(.01f);
+
+ cannyEdgeDetectionfilter->Update();
+ m_EdgeImage = cannyEdgeDetectionfilter->GetOutput();
+
+
+ // set minCosts
+ minCosts = 0.0; // The lower, the more thouroughly! 0 = dijkstra. If estimate costs are lower than actual costs everything is fine. If estimation is higher than actual costs, you might not get the shortest but a different path.
+
+ m_Initialized = true;
+ }
+
+ // check start/end point value
+ startValue= this->m_Image->GetPixel(this->m_StartIndex);
+ endValue= this->m_Image->GetPixel(this->m_EndIndex);
+ }
+
+
+
+ template
+ double ShortestPathCostFunctionLiveWire::SigmoidFunction(double I, double max, double min, double alpha, double beta)
+ {
+ // Using the SIgmoid formula from ITK Software Guide 6.3.2 Non Linear Mappings
+ double Exponent = -1 * ((I - beta) / alpha);
+ double Factor = 1 / (1 + exp(Exponent));
+ double newI = (max - min) * Factor + min;
+
+ return newI;
+ }
+
+
+ template
+ double ShortestPathCostFunctionLiveWire::Gaussian(double x, double xOfGaussian, double yOfGaussian)
+ {
+ return yOfGaussian * exp( -0.5 * pow( (x - xOfGaussian), 2) );
+ }
+
+} // end namespace itk
+
+#endif // __itkShortestPathCostFunctionLiveWire_txx
diff --git a/Modules/GraphAlgorithms/itkShortestPathImageFilter.h b/Modules/GraphAlgorithms/itkShortestPathImageFilter.h
index c93a0204a3..1a3216a6dc 100644
--- a/Modules/GraphAlgorithms/itkShortestPathImageFilter.h
+++ b/Modules/GraphAlgorithms/itkShortestPathImageFilter.h
@@ -1,225 +1,232 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef __itkShortestPathImageFilter_h
#define __itkShortestPathImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkShortestPathCostFunction.h"
#include "itkShortestPathNode.h"
#include
#include
// ------- INFORMATION ----------
/// SET FUNCTIONS
//void SetInput( ItkImage ) // Compulsory
//void SetStartIndex (const IndexType & StartIndex); // Compulsory
//void SetEndIndex(const IndexType & EndIndex); // Compulsory
//void SetFullNeighborsMode(bool) // Optional (default=false), if false N4, if true N26
//void SetActivateTimeOut(bool) // Optional (default=false), for debug issues: after 30s algorithms terminates. You can have a look at the VectorOrderImage to see how far it came
//void SetMakeOutputImage(bool) // Optional (default=true), Generate an outputimage of the path. You can also get the path directoy with GetVectorPath()
//void SetCalcAllDistances(bool) // Optional (default=false), Calculate Distances over the whole image. CAREFUL, algorithm time extends a lot. Necessary for GetDistanceImage
//void SetStoreVectorOrder(bool) // Optional (default=false), Stores in which order the pixels were checked. Necessary for GetVectorOrderImage
//void AddEndIndex(const IndexType & EndIndex) //Optional. By calling this function you can add several endpoints! The algorithm will look for several shortest Pathes. From Start to all Endpoints.
//
/// GET FUNCTIONS
//std::vector< itk::Index<3> > GetVectorPath(); // returns the shortest path as vector
//std::vector< std::vector< itk::Index<3> > GetMultipleVectorPathe(); // returns a vector of shortest Pathes (which are vectors of points)
//GetDistanceImage // Returns the distance image
//GetVectorOrderIMage // Returns the Vector Order image
//
// EXAMPLE USE
// pleae see qmitkmitralvalvesegmentation4dtee bundle
namespace itk
{
template
class ShortestPathImageFilter :
public ImageToImageFilter
{
public:
//Standard Typedefs
typedef ShortestPathImageFilter Self;
typedef ImageToImageFilter Superclass;
typedef SmartPointer Pointer;
typedef SmartPointer ConstPointer;
// Typdefs for metric
typedef ShortestPathCostFunction< TInputImageType > CostFunctionType;
typedef typename CostFunctionType::Pointer CostFunctionTypePointer;
// More typdefs for convenience
typedef TInputImageType InputImageType;
typedef typename TInputImageType::Pointer InputImagePointer;
typedef typename TInputImageType::PixelType InputImagePixelType;
typedef typename TInputImageType::SizeType InputImageSizeType;
typedef typename TInputImageType::IndexType IndexType;
typedef typename itk::ImageRegionIteratorWithIndex< InputImageType > InputImageIteratorType;
typedef TOutputImageType OutputImageType;
typedef typename TOutputImageType::Pointer OutputImagePointer;
typedef typename TOutputImageType::PixelType OutputImagePixelType;
typedef typename TOutputImageType::IndexType OutputImageIndexType;
typedef ImageRegionIteratorWithIndex< OutputImageType > OutputImageIteratorType;
typedef itk::ShapedNeighborhoodIterator< TInputImageType > itkShapedNeighborhoodIteratorType;
// New Macro for smartpointer instantiation
itkNewMacro(Self)
// Run-time type information
itkTypeMacro(ShortestPathImageFilter, ImageToImageFilter)
// Display
void PrintSelf( std::ostream& os, Indent indent ) const;
// Compare function for A_STAR
struct CompareNodeStar
{
bool operator()(ShortestPathNode *a, ShortestPathNode *b)
{
return (a->distAndEst > b->distAndEst);
}
};
// \brief Set Starpoint for ShortestPath Calculation
void SetStartIndex (const IndexType & StartIndex);
// \brief Adds Endpoint for multiple ShortestPath Calculation
void AddEndIndex (const IndexType & index);
// \brief Set Endpoint for ShortestPath Calculation
void SetEndIndex(const IndexType & EndIndex);
// \brief Set FullNeighborsMode. false = no diagonal neighbors, in 2D this means N4 Neigborhood. true = would be N8 in 2D
itkSetMacro (FullNeighborsMode, bool);
itkGetMacro (FullNeighborsMode, bool);
// \brief (default=true), Produce output image, which shows the shortest path. But you can also get the shortest Path directly as vector with the function GetVectorPath
itkSetMacro (MakeOutputImage, bool);
itkGetMacro (MakeOutputImage, bool);
// \brief (default=false), Store an Vector of Order, so you can call getVectorOrderImage after update
itkSetMacro (StoreVectorOrder, bool);
itkGetMacro (StoreVectorOrder, bool);
// \brief (default=false), // Calculate all Distances to all pixels, so you can call getDistanceImage after update (warning algo will take a long time)
itkSetMacro (CalcAllDistances, bool);
itkGetMacro (CalcAllDistances, bool);
// \brief (default=false), for debug issues: after 30s algorithms terminates. You can have a look at the VectorOrderImage to see how far it came
itkSetMacro (ActivateTimeOut, bool);
itkGetMacro (ActivateTimeOut, bool);
// \brief returns shortest Path as vector
- std::vector< itk::Index<3> > GetVectorPath();
+ std::vector< IndexType > GetVectorPath();
// \brief returns Multiple shortest Paths. You can call this function, when u performed a multiple shortest path search (one start, several ends)
- std::vector< std::vector< itk::Index<3> > > GetMultipleVectorPaths();
+ std::vector< std::vector< IndexType > > GetMultipleVectorPaths();
// \brief returns the vector order image. It shows in which order the pixels were checked. good for debugging. Be sure to have m_StoreVectorOrder=true
OutputImagePointer GetVectorOrderImage();
// \brief returns the distance image. It shows the distances from the startpoint to all other pixels. Be sure to have m_CalcAllDistances=true
OutputImagePointer GetDistanceImage();
+
+ // \brief Fill m_VectorPath
+ void MakeShortestPathVector();
+
// \brief cleans up the filter
void CleanUp();
itkSetObjectMacro( CostFunction, CostFunctionType ); // itkSetObjectMacro = set function that uses pointer as parameter
itkGetObjectMacro( CostFunction, CostFunctionType );
protected:
std::vector< IndexType > m_endPoints; // if you fill this vector, the algo will not rest until all endPoints have been reached
std::vector< IndexType > m_endPointsClosed;
ShortestPathNode* m_Nodes; // main list that contains all nodes
NodeNumType m_Graph_NumberOfNodes;
NodeNumType m_Graph_StartNode;
NodeNumType m_Graph_EndNode;
int m_ImageDimensions;
bool m_Graph_fullNeighbors;
std::vector m_Graph_DiscoveredNodeList;
ShortestPathImageFilter(Self&); // intentionally not implemented
void operator=(const Self&); // intentionally not implemented
const static int BACKGROUND = 0;
const static int FOREGROUND = 255;
bool m_FullNeighborsMode;
bool m_MakeOutputImage;
bool m_StoreVectorOrder; // Store an Vector of Order, so you can call getVectorOrderImage after update
bool m_CalcAllDistances; // Calculate all Distances, so you can call getDistanceImage after update (warning algo will take a long time)
bool multipleEndPoints;
bool m_ActivateTimeOut; // if true, then i search max. 30 secs. then abort
+
+ bool m_Initialized;
+
+
CostFunctionTypePointer m_CostFunction;
IndexType m_StartIndex, m_EndIndex;
- std::vector< itk::Index<3> > m_VectorPath;
- std::vector< std::vector< itk::Index<3> > > m_MultipleVectorPaths;
+ std::vector< IndexType > m_VectorPath;
+ std::vector< std::vector< IndexType > > m_MultipleVectorPaths;
std::vector< NodeNumType > m_VectorOrder;
ShortestPathImageFilter();
- // \brief Fill m_VectorPath
- void MakeShortestPathVector();
+ ~ShortestPathImageFilter();
// \brief Create all the outputs
void MakeOutputs();
// \brief Generate Data
void GenerateData();
// \brief gets the estimate costs from pixel a to target.
double getEstimatedCostsToTarget(const IndexType & a);
typename InputImageType::Pointer m_magnitudeImage;
// \brief Convert a indexnumber of a node in m_Nodes to image coordinates
typename TInputImageType::IndexType NodeToCoord(NodeNumType);
// \brief Convert image coordinate to a indexnumber of a node in m_Nodes
unsigned int CoordToNode(IndexType);
// \brief Returns the neighbors of a node
std::vector GetNeighbors(NodeNumType nodeNum, bool FullNeighbors);
// \brief Check if coords are in bounds of image
bool CoordIsInBounds(IndexType);
// \brief Initializes the graph
void InitGraph();
// \brief Start ShortestPathSearch
void StartShortestPathSearch();
};
} // end of namespace itk
#include "itkShortestPathImageFilter.txx"
#endif
diff --git a/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx b/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx
index a255040bfb..933130fe43 100644
--- a/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx
+++ b/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx
@@ -1,939 +1,951 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef __itkShortestPathImageFilter_txx
#define __itkShortestPathImageFilter_txx
#include "time.h"
#include "mitkMemoryUtilities.h"
#include
#include
#include
namespace itk
{
// Constructor (initialize standard values)
template
ShortestPathImageFilter
::ShortestPathImageFilter() :
m_FullNeighborsMode(false),
m_MakeOutputImage(true),
m_StoreVectorOrder(false),
m_CalcAllDistances(false),
m_ActivateTimeOut(false),
- multipleEndPoints(false),
- m_Nodes(0),
- m_Graph_NumberOfNodes(0)
+ multipleEndPoints(false),
+ m_Initialized(false),
+ m_Nodes(0),
+ m_Graph_NumberOfNodes(0)
{
- m_endPoints.clear();
- m_endPointsClosed.clear();
+ m_endPoints.clear();
+ m_endPointsClosed.clear();
- if (m_MakeOutputImage)
- {
- this->SetNumberOfRequiredOutputs(1);
- this->SetNthOutput( 0, OutputImageType::New() );
- }
+ if (m_MakeOutputImage)
+ {
+ this->SetNumberOfRequiredOutputs(1);
+ this->SetNthOutput( 0, OutputImageType::New() );
+ }
}
+
+
+ template
+ ShortestPathImageFilter
+ ::~ShortestPathImageFilter()
+ {
+ delete [] m_Nodes;
+ }
+
+
+
+
template
inline typename ShortestPathImageFilter::IndexType
ShortestPathImageFilter
::NodeToCoord (NodeNumType node)
{
const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize();
int dim = InputImageType::ImageDimension;
IndexType coord;
if (dim == 2)
{
coord[1] = node / size[0];
coord[0] = node % size[0];
if ((coord[0] >= size[0]) || (coord[1] >= size[1]))
{
coord[0] = 0;
coord[1] = 0;
}
}
if (dim == 3)
{
coord[2] = node / (size[0]*size[1]);
coord[1] = (node % (size[0]*size[1])) / size[0];
coord[0] = (node % (size[0]*size[1])) % size[0];
if ((coord[0] >= size[0]) || (coord[1] >= size[1]) || (coord[2] >= size[2]))
{
coord[0] = 0;
coord[1] = 0;
coord[2] = 0;
}
}
return coord;
}
template
inline typename itk::NodeNumType
ShortestPathImageFilter::
CoordToNode (IndexType coord)
{
const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize();
int dim = InputImageType::ImageDimension;
NodeNumType node = 0;
if (dim == 2)
{
node = (coord[1]*size[0]) + coord[0];
}
if (dim == 3)
{
node = (coord[2]*size[0]*size[1]) + (coord[1]*size[0]) + coord[0];
}
if ((m_Graph_NumberOfNodes > 0) && (node >= m_Graph_NumberOfNodes))
- {
- MITK_INFO << "WARNING! Coordinates outside image!";
- MITK_INFO << "Coords = " << coord ;
- MITK_INFO << "ImageDim = " << dim ;
- MITK_INFO << "RequestedRegionSize = " << size ;
+ {
+ /*MITK_INFO << "WARNING! Coordinates outside image!";
+ MITK_INFO << "Coords = " << coord ;
+ MITK_INFO << "ImageDim = " << dim ;
+ MITK_INFO << "RequestedRegionSize = " << size ;*/
node = 0;
- }
+ }
return node;
}
template
inline bool
ShortestPathImageFilter::
CoordIsInBounds (IndexType coord)
{
const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize();
int dim = InputImageType::ImageDimension;
if (dim == 2)
{
if ((coord[0] >= 0)
&& (coord[0] < size[0])
&& (coord[1] >= 0 )
&& (coord[1] < size[1] ))
{
return true;
}
}
if (dim == 3)
{
if ((coord[0] >= 0)
&& (coord[0] < size[0])
&& (coord[1] >= 0 )
&& (coord[1] < size[1] )
&& (coord[2] >= 0 )
&& (coord[2] < size[2] ))
{
return true;
}
}
return false;
}
template
inline std::vector< ShortestPathNode* >
ShortestPathImageFilter::
GetNeighbors (unsigned int nodeNum, bool FullNeighbors)
{
// returns a vector of nodepointers.. these nodes are the neighbors
int dim = InputImageType::ImageDimension;
IndexType Coord = NodeToCoord(nodeNum);
IndexType NeighborCoord;
std::vector nodeList;
int neighborDistance = 1; //if i increase that, i might not hit the endnote
- // maybe use itkNeighborhoodIterator here, might be faster
+ // maybe use itkNeighborhoodIterator here, might be faster
if ( dim == 2)
{
// N4
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
if (FullNeighbors)
{
// N8
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
}
}
if ( dim == 3)
{
// N6
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
if (FullNeighbors)
{
// N26
// Middle Slice
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2];
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
// BackSlice (Diagonal)
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
//BackSlice (Non-Diag)
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]-neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
// FrontSlice (Diagonal)
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
//FrontSlice(Non-Diag)
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]-neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]+neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0];
NeighborCoord[1] = Coord[1]+neighborDistance;
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
NeighborCoord[0] = Coord[0]-neighborDistance;
NeighborCoord[1] = Coord[1];
NeighborCoord[2] = Coord[2]+neighborDistance;
if (CoordIsInBounds(NeighborCoord))
nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]);
}
}
return nodeList;
}
template
void ShortestPathImageFilter::
SetStartIndex (const typename TInputImageType::IndexType &StartIndex)
{
for (unsigned int i=0;i
void ShortestPathImageFilter::
SetEndIndex (const typename TInputImageType::IndexType &EndIndex)
{
for (unsigned int i=0;i
+ template
void ShortestPathImageFilter::
- AddEndIndex (const typename TInputImageType::IndexType &index)
+ AddEndIndex (const typename TInputImageType::IndexType &index)
{
- // ONLY FOR MULTIPLE END POINTS SEARCH
- IndexType newEndIndex;
+ // ONLY FOR MULTIPLE END POINTS SEARCH
+ IndexType newEndIndex;
for (unsigned int i=0;i
inline double ShortestPathImageFilter::
getEstimatedCostsToTarget (const typename TInputImageType::IndexType &a)
{
// Returns the minimal possible costs for a path from "a" to targetnode.
- itk::Vector v;
+ itk::Vector v;
v[0] = m_EndIndex[0]-a[0];
v[1] = m_EndIndex[1]-a[1];
v[2] = m_EndIndex[2]-a[2];
- return m_CostFunction->GetMinCost() * v.GetNorm();
+ return m_CostFunction->GetMinCost() * v.GetNorm();
}
template
void
ShortestPathImageFilter::
InitGraph()
{
- // Clean up previous stuff
- CleanUp();
+ if(!m_Initialized)
+ {
+ // Clean up previous stuff
+ CleanUp();
- // initalize cost function
- m_CostFunction->Initialize();
+ // Calc Number of nodes
+ m_ImageDimensions = TInputImageType::ImageDimension;
+ const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize();
+ m_Graph_NumberOfNodes = 1;
+ for (NodeNumType i=0; iGetInput()->GetRequestedRegion().GetSize();
- m_Graph_NumberOfNodes = 1;
- for (NodeNumType i=0; iInitialize();
}
template
void
ShortestPathImageFilter::
StartShortestPathSearch()
{
- // Setup Timer
+ // Setup Timer
clock_t startAll = clock();
clock_t stopAll = clock();
- // init variables
+ // init variables
double durationAll = 0;
bool timeout = false;
bool makeNewHeapNecessary = false;
- NodeNumType mainNodeListIndex = 0;
+ NodeNumType mainNodeListIndex = 0;
DistanceType curNodeDistance = 0;
DistanceType curNodeDistAndEst = 0;
- NodeNumType numberOfNodesChecked = 0;
+ NodeNumType numberOfNodesChecked = 0;
- // Create Multimap (tree structure for fast searching)
- std::multimap myMap;
- std::pair< std::multimap::iterator, std::multimap::iterator> ret;
- std::multimap::iterator it;
+ // Create Multimap (tree structure for fast searching)
+ std::multimap myMap;
+ std::pair< std::multimap::iterator, std::multimap::iterator> ret;
+ std::multimap::iterator it;
- // At first, only startNote is discovered.
- myMap.insert( std::pair (m_Nodes[m_Graph_StartNode].distAndEst, &m_Nodes[m_Graph_StartNode]) );
+ // At first, only startNote is discovered.
+ myMap.insert( std::pair (m_Nodes[m_Graph_StartNode].distAndEst, &m_Nodes[m_Graph_StartNode]) );
// While there are discovered Nodes, pick the one with lowest distance,
- // update its neighbors and eventually delete it from the discovered Nodes list.
+ // update its neighbors and eventually delete it from the discovered Nodes list.
while(!myMap.empty())
{
- numberOfNodesChecked++;
- if ( (numberOfNodesChecked % (m_Graph_NumberOfNodes/100)) == 0)
- {
- MITK_INFO << "Checked " << ( numberOfNodesChecked / (m_Graph_NumberOfNodes/100) ) << "% List: " << myMap.size() << "\n";
- }
-
- // Get element with lowest score
- mainNodeListIndex = myMap.begin()->second->mainListIndex;
+ numberOfNodesChecked++;
+ /*if ( (numberOfNodesChecked % (m_Graph_NumberOfNodes/100)) == 0)
+ {
+ MITK_INFO << "Checked " << ( numberOfNodesChecked / (m_Graph_NumberOfNodes/100) ) << "% List: " << myMap.size() << "\n";
+ }*/
+
+ // Get element with lowest score
+ mainNodeListIndex = myMap.begin()->second->mainListIndex;
curNodeDistAndEst = myMap.begin()->second->distAndEst;
curNodeDistance = myMap.begin()->second->distance;
- myMap.begin()->second->closed = true; // close it
+ myMap.begin()->second->closed = true; // close it
- // Debug:
- //MITK_INFO << "INFO: size " << myMap.size();
- /*
- for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
- }
- */
+ // Debug:
+ //MITK_INFO << "INFO: size " << myMap.size();
+ /*
+ for (it = myMap.begin(); it != myMap.end(); ++it)
+ {
+ MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
+ }
+ */
// Kicks out element with lowest score
- myMap.erase( myMap.begin() );
+ myMap.erase( myMap.begin() );
- // if wanted, store vector order
+ // if wanted, store vector order
if (m_StoreVectorOrder)
- {
+ {
m_VectorOrder.push_back(mainNodeListIndex);
- }
+ }
// Check neighbors
std::vector neighborNodes = GetNeighbors(mainNodeListIndex, m_Graph_fullNeighbors);
for (NodeNumType i=0; iclosed)
- continue; // this nodes is already closed, go to next neighbor
+ if (neighborNodes[i]->closed)
+ continue; // this nodes is already closed, go to next neighbor
IndexType coordCurNode = NodeToCoord(mainNodeListIndex);
IndexType coordNeighborNode = NodeToCoord(neighborNodes[i]->mainListIndex);
// calculate the new Distance to the current neighbor
double newDistance = curNodeDistance
+ (m_CostFunction->GetCost(coordCurNode, coordNeighborNode));
// if it is shorter than any yet known path to this neighbor, than the current path is better. Save that!
if ((newDistance < neighborNodes[i]->distance) || (neighborNodes[i]->distance == -1) )
{
// if that neighbornode is not in discoverednodeList yet, Push it there and update
if (neighborNodes[i]->distance == -1)
{
neighborNodes[i]->distance = newDistance;
neighborNodes[i]->distAndEst = newDistance + getEstimatedCostsToTarget(coordNeighborNode);
neighborNodes[i]->prevNode = mainNodeListIndex;
- myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) );
- /*
- MITK_INFO << "Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex;
- MITK_INFO << "INFO: size " << myMap.size();
- for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
- }
- */
+ myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) );
+ /*
+ MITK_INFO << "Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex;
+ MITK_INFO << "INFO: size " << myMap.size();
+ for (it = myMap.begin(); it != myMap.end(); ++it)
+ {
+ MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
+ }
+ */
}
// or if is already in discoverednodelist, update
else
{
- /*
- it = myMap.find(neighborNodes[i]->distAndEst);
- if (it == myMap.end() )
- {
- MITK_INFO << "Nothing!";
- // look further
- for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- if ((*it).second->mainListIndex == lookForId)
- {
- MITK_INFO << "But it is there!!!";
- MITK_INFO << "Searched for: " << lookFor << " but had: " << (*it).second->distAndEst;
- }
-
- }
- }
- */
+ /*
+ it = myMap.find(neighborNodes[i]->distAndEst);
+ if (it == myMap.end() )
+ {
+ MITK_INFO << "Nothing!";
+ // look further
+ for (it = myMap.begin(); it != myMap.end(); ++it)
+ {
+ if ((*it).second->mainListIndex == lookForId)
+ {
+ MITK_INFO << "But it is there!!!";
+ MITK_INFO << "Searched for: " << lookFor << " but had: " << (*it).second->distAndEst;
+ }
- // 1st : find and delete old element
- bool found = false;
- double lookFor = neighborNodes[i]->distAndEst;
- unsigned int lookForId = neighborNodes[i]->mainListIndex;
- ret = myMap.equal_range(neighborNodes[i]->distAndEst);
+ }
+ }
+ */
- if ((ret.first == ret.second))
- {
- MITK_INFO << "No exact match!"; // if this happens, you are screwed
- /*
- MITK_INFO << "Was looking for: " << lookFor << " ID: " << lookForId;
- if (ret.first != myMap.end() )
- {
- it = ret.first;
- MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
- ++it;
- MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
- --it;
- --it;
- MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
- }
-
- // look if that ID is found in the map
- for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- if ((*it).second->mainListIndex == lookForId)
- {
- MITK_INFO << "But it is there!!!";
- MITK_INFO << "Searched dist: " << lookFor << " found dist: " << (*it).second->distAndEst;
- }
-
- }
- */
- }
- else
- {
- for (it=ret.first; it!=ret.second; ++it)
- {
- if (it->second->mainListIndex == neighborNodes[i]->mainListIndex)
- {
- found = true;
- myMap.erase(it);
- /*
- MITK_INFO << "INFO: size " << myMap.size();
- MITK_INFO << "Erase: " << it->first << "|" << it->second->mainListIndex;
-
- MITK_INFO << "INFO: size " << myMap.size();
- for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
- }
- */
- break;
- }
- }
- }
+ // 1st : find and delete old element
+ bool found = false;
+ double lookFor = neighborNodes[i]->distAndEst;
+ unsigned int lookForId = neighborNodes[i]->mainListIndex;
+ ret = myMap.equal_range(neighborNodes[i]->distAndEst);
+ if ((ret.first == ret.second))
+ {
+ /*+++++++++++++ no exact match +++++++++++++*/
+ //MITK_INFO << "No exact match!"; // if this happens, you are screwed
+ /*
+ MITK_INFO << "Was looking for: " << lookFor << " ID: " << lookForId;
+ if (ret.first != myMap.end() )
+ {
+ it = ret.first;
+ MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
+ ++it;
+ MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
+ --it;
+ --it;
+ MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex;
+ }
+
+ // look if that ID is found in the map
+ for (it = myMap.begin(); it != myMap.end(); ++it)
+ {
+ if ((*it).second->mainListIndex == lookForId)
+ {
+ MITK_INFO << "But it is there!!!";
+ MITK_INFO << "Searched dist: " << lookFor << " found dist: " << (*it).second->distAndEst;
+ }
+
+ }
+ */
+ }
+ else
+ {
+ for (it=ret.first; it!=ret.second; ++it)
+ {
+ if (it->second->mainListIndex == neighborNodes[i]->mainListIndex)
+ {
+ found = true;
+ myMap.erase(it);
+ /*
+ MITK_INFO << "INFO: size " << myMap.size();
+ MITK_INFO << "Erase: " << it->first << "|" << it->second->mainListIndex;
- if (!found)
+ MITK_INFO << "INFO: size " << myMap.size();
+ for (it = myMap.begin(); it != myMap.end(); ++it)
{
- MITK_INFO << "Could not find it! :(";
- continue;
+ MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
}
+ */
+ break;
+ }
+ }
+ }
- // 2nd: update and insert new element
- neighborNodes[i]->distance = newDistance;
+ if (!found)
+ {
+ //MITK_INFO << "Could not find it! :(";
+ continue;
+ }
+
+ // 2nd: update and insert new element
+ neighborNodes[i]->distance = newDistance;
neighborNodes[i]->distAndEst = newDistance + getEstimatedCostsToTarget(coordNeighborNode);
neighborNodes[i]->prevNode = mainNodeListIndex;
- //myMap.insert( std::pair (neighborNodes[i]->distAndEst, neighborNodes[i]));
- myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) );
+ //myMap.insert( std::pair (neighborNodes[i]->distAndEst, neighborNodes[i]));
+ myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) );
- //MITK_INFO << "Re-Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex;
- //MITK_INFO << "INFO: size " << myMap.size();
- /*for (it = myMap.begin(); it != myMap.end(); ++it)
- {
- MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
- }*/
+ //MITK_INFO << "Re-Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex;
+ //MITK_INFO << "INFO: size " << myMap.size();
+ /*for (it = myMap.begin(); it != myMap.end(); ++it)
+ {
+ MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex;
+ }*/
}
}
}
// finished with checking all neighbors.
// Check Timeout, if activated
if (m_ActivateTimeOut)
{
stopAll = clock();
durationAll = (double)(stopAll - startAll) / CLOCKS_PER_SEC;
if (durationAll >= 30)
{
- MITK_INFO << "TIMEOUT!! Search took over 30 seconds";
+ //MITK_INFO << "TIMEOUT!! Search took over 30 seconds";
timeout = true ;
}
}
- // Check end criteria:
- // For multiple points
- if ( multipleEndPoints )
- {
- // super slow, make it faster
- for (int i=0 ;i
void
ShortestPathImageFilter::
MakeOutputs()
{
- MITK_INFO << "Make Output";
+ //MITK_INFO << "Make Output";
- if (m_MakeOutputImage)
+ if (m_MakeOutputImage)
+ {
+ OutputImagePointer output0 = this->GetOutput(0);
+ output0->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
+ output0->Allocate();
+ OutputImageIteratorType shortestPathImageIt (output0, output0->GetRequestedRegion());
+ // Create ShortestPathImage (Output 0)
+ for (shortestPathImageIt.GoToBegin(); !shortestPathImageIt.IsAtEnd(); ++shortestPathImageIt)
{
- OutputImagePointer output0 = this->GetOutput(0);
- output0->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
- output0->Allocate();
- OutputImageIteratorType shortestPathImageIt (output0, output0->GetRequestedRegion());
- // Create ShortestPathImage (Output 0)
- for (shortestPathImageIt.GoToBegin(); !shortestPathImageIt.IsAtEnd(); ++shortestPathImageIt)
- {
- // First intialize with background color
- shortestPathImageIt.Set( BACKGROUND ) ;
- }
-
- if (!multipleEndPoints) // Only one path was calculated
- {
- for (int i=0; i< m_VectorPath.size(); i++)
+ // First intialize with background color
+ shortestPathImageIt.Set( BACKGROUND ) ;
+ }
+
+ if (!multipleEndPoints) // Only one path was calculated
+ {
+ for (int i=0; i< m_VectorPath.size(); i++)
+ {
+ shortestPathImageIt.SetIndex( m_VectorPath[i] );
+ shortestPathImageIt.Set( FOREGROUND ) ;
+ }
+ }
+ else // multiple pathes has been calculated, draw all
+ {
+ for (int i =0; i
- typename ShortestPathImageFilter::OutputImagePointer
+ typename ShortestPathImageFilter::OutputImagePointer
ShortestPathImageFilter::
GetVectorOrderImage()
{
- // Create Vector Order Image
- // Return it
+ // Create Vector Order Image
+ // Return it
- OutputImagePointer image = OutputImageType::New();
- image->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
- image->Allocate();
- OutputImageIteratorType vectorOrderImageIt (image, image->GetRequestedRegion());
+ OutputImagePointer image = OutputImageType::New();
+ image->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
+ image->Allocate();
+ OutputImageIteratorType vectorOrderImageIt (image, image->GetRequestedRegion());
- MITK_INFO << "GetVectorOrderImage";
- for (vectorOrderImageIt.GoToBegin(); !vectorOrderImageIt.IsAtEnd(); ++vectorOrderImageIt)
- {
- // First intialize with background color
- vectorOrderImageIt.Value() = BACKGROUND ;
- }
- for (int i=0; i< m_VectorOrder.size(); i++)
- {
- vectorOrderImageIt.SetIndex(NodeToCoord(m_VectorOrder[i]) );
- vectorOrderImageIt.Set( BACKGROUND+i) ;
- }
+ //MITK_INFO << "GetVectorOrderImage";
+ for (vectorOrderImageIt.GoToBegin(); !vectorOrderImageIt.IsAtEnd(); ++vectorOrderImageIt)
+ {
+ // First intialize with background color
+ vectorOrderImageIt.Value() = BACKGROUND ;
+ }
+ for (int i=0; i< m_VectorOrder.size(); i++)
+ {
+ vectorOrderImageIt.SetIndex(NodeToCoord(m_VectorOrder[i]) );
+ vectorOrderImageIt.Set( BACKGROUND+i) ;
+ }
return image;
}
template
typename ShortestPathImageFilter::OutputImagePointer
ShortestPathImageFilter::
GetDistanceImage()
{
- // Create Distance Image
- // Return it
-
- OutputImagePointer image = OutputImageType::New();
- image->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
- image->Allocate();;
- OutputImageIteratorType distanceImageIt (image, image->GetRequestedRegion());
- // Create Distance Image (Output 1)
- NodeNumType myNodeNum;
- for (distanceImageIt.GoToBegin(); !distanceImageIt.IsAtEnd(); ++distanceImageIt)
- {
- IndexType index = distanceImageIt.GetIndex();
- myNodeNum = CoordToNode(index);
- double newVal = m_Nodes[myNodeNum].distance;
- distanceImageIt.Set(newVal);
- }
- }
+ // Create Distance Image
+ // Return it
+
+ OutputImagePointer image = OutputImageType::New();
+ image->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
+ image->Allocate();;
+ OutputImageIteratorType distanceImageIt (image, image->GetRequestedRegion());
+ // Create Distance Image (Output 1)
+ NodeNumType myNodeNum;
+ for (distanceImageIt.GoToBegin(); !distanceImageIt.IsAtEnd(); ++distanceImageIt)
+ {
+ IndexType index = distanceImageIt.GetIndex();
+ myNodeNum = CoordToNode(index);
+ double newVal = m_Nodes[myNodeNum].distance;
+ distanceImageIt.Set(newVal);
+ }
+ }
template
- std::vector< itk::Index<3> >
+ std::vector< typename ShortestPathImageFilter::IndexType >
ShortestPathImageFilter::
GetVectorPath()
{
return m_VectorPath;
}
- template
- std::vector< std::vector< itk::Index<3> > >
+ template
+ std::vector< std::vector< typename ShortestPathImageFilter::IndexType > >
ShortestPathImageFilter::
GetMultipleVectorPaths()
{
return m_MultipleVectorPaths;
}
template
void
ShortestPathImageFilter::
MakeShortestPathVector()
{
- MITK_INFO << "Make ShortestPath Vec";
-
- // single end point
- if ( !multipleEndPoints )
- {
- // fill m_VectorPath with the Shortest Path
- m_VectorPath.clear();
+ //MITK_INFO << "Make ShortestPath Vec";
- // Go backwards from endnote to startnode
- NodeNumType prevNode = m_Graph_EndNode;
- while (prevNode != -1)
- {
- m_VectorPath.push_back( NodeToCoord(prevNode) );
+ // single end point
+ if ( !multipleEndPoints )
+ {
+ // fill m_VectorPath with the Shortest Path
+ m_VectorPath.clear();
- if (prevNode == m_Graph_StartNode)
- break;
+ // Go backwards from endnote to startnode
+ NodeNumType prevNode = m_Graph_EndNode;
+ while (prevNode != -1)
+ {
+ m_VectorPath.push_back( NodeToCoord(prevNode) );
- prevNode = m_Nodes[prevNode].prevNode;
- }
+ if (prevNode == m_Graph_StartNode)
+ break;
- // reverse it
- std::reverse(m_VectorPath.begin(), m_VectorPath.end() );
+ prevNode = m_Nodes[prevNode].prevNode;
}
- // Multiple end end points and pathes
- else
+
+ // reverse it
+ std::reverse(m_VectorPath.begin(), m_VectorPath.end() );
+ }
+ // Multiple end end points and pathes
+ else
+ {
+ for (int i=0; i
void
ShortestPathImageFilter::
CleanUp()
{
-
- m_VectorPath.clear();
- //TODO: if multiple Path, clear all multiple Paths
- /*
- for (NodeNumType i=0; i
void
ShortestPathImageFilter::
GenerateData()
{
// Build Graph
InitGraph();
// Calc Shortest Parth
StartShortestPathSearch();
// Fill Shortest Path
MakeShortestPathVector();
// Make Outputs
MakeOutputs();
}
template
void
ShortestPathImageFilter::
PrintSelf( std::ostream& os, Indent indent ) const
{
Superclass::PrintSelf(os,indent);
}
} /* end namespace itk */
#endif
diff --git a/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.cpp b/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.cpp
new file mode 100644
index 0000000000..69094bfa56
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.cpp
@@ -0,0 +1,416 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include "mitkImageLiveWireContourModelFilter.h"
+
+
+#include
+#include
+#include
+
+
+mitk::ImageLiveWireContourModelFilter::ImageLiveWireContourModelFilter()
+{
+ OutputType::Pointer output = dynamic_cast ( this->MakeOutput( 0 ).GetPointer() );
+ this->SetNumberOfRequiredInputs(1);
+ this->SetNumberOfOutputs( 1 );
+ this->SetNthOutput(0, output.GetPointer());
+ m_CostFunction = ImageLiveWireContourModelFilter::CostFunctionType::New();
+ m_ShortestPathFilter = ShortestPathImageFilterType::New();
+ m_ShortestPathFilter->SetCostFunction(m_CostFunction);
+ m_UseDynamicCostMap = false;
+ m_ImageModified = false;
+ m_Timestep = 0;
+}
+
+mitk::ImageLiveWireContourModelFilter::~ImageLiveWireContourModelFilter()
+{
+
+}
+
+
+mitk::ImageLiveWireContourModelFilter::OutputType* mitk::ImageLiveWireContourModelFilter::GetOutput()
+{
+ return Superclass::GetOutput();
+}
+
+void mitk::ImageLiveWireContourModelFilter::SetInput ( const mitk::ImageLiveWireContourModelFilter::InputType* input )
+{
+ this->SetInput( 0, input );
+}
+
+void mitk::ImageLiveWireContourModelFilter::SetInput ( unsigned int idx, const mitk::ImageLiveWireContourModelFilter::InputType* input )
+{
+ if ( idx + 1 > this->GetNumberOfInputs() )
+ {
+ this->SetNumberOfRequiredInputs(idx + 1);
+ }
+ if ( input != static_cast ( this->ProcessObject::GetInput ( idx ) ) )
+ {
+ this->ProcessObject::SetNthInput ( idx, const_cast ( input ) );
+ this->Modified();
+ this->m_ImageModified = true;
+ m_ShortestPathFilter = ShortestPathImageFilterType::New();
+ m_ShortestPathFilter->SetCostFunction(m_CostFunction);
+ }
+}
+
+
+
+const mitk::ImageLiveWireContourModelFilter::InputType* mitk::ImageLiveWireContourModelFilter::GetInput( void )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(0));
+}
+
+
+const mitk::ImageLiveWireContourModelFilter::InputType* mitk::ImageLiveWireContourModelFilter::GetInput( unsigned int idx )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(idx));
+}
+
+
+void mitk::ImageLiveWireContourModelFilter::GenerateData()
+{
+ mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+
+ if(!input)
+ {
+ MITK_ERROR << "No input available.";
+ itkExceptionMacro("mitk::ImageToLiveWireContourFilter: No input available. Please set the input!");
+ return;
+ }
+
+ if( input->GetDimension() != 2 )
+ {
+ MITK_ERROR << "Filter is only working on 2D images.";
+ itkExceptionMacro("mitk::ImageToLiveWireContourFilter: Filter is only working on 2D images.. Please make sure that the input is 2D!");
+ return;
+ }
+
+
+ input->GetGeometry()->WorldToIndex(m_StartPoint, m_StartPointInIndex);
+ input->GetGeometry()->WorldToIndex(m_EndPoint, m_EndPointInIndex);
+
+ //only start calculating if both indices are inside image geometry
+ if( input->GetGeometry()->IsIndexInside(this->m_StartPointInIndex) && input->GetGeometry()->IsIndexInside(this->m_EndPointInIndex) )
+ {
+ AccessFixedDimensionByItk(input, ItkProcessImage, 2);
+ m_ImageModified = false;
+ }
+}
+
+
+
+template
+void mitk::ImageLiveWireContourModelFilter::ItkProcessImage (itk::Image* inputImage)
+{
+ typedef itk::Image< TPixel, VImageDimension > InputImageType;
+ typedef typename InputImageType::IndexType IndexType;
+
+
+ /* compute the requested region for itk filters */
+
+ IndexType startPoint, endPoint;
+
+ startPoint[0] = m_StartPointInIndex[0];
+ startPoint[1] = m_StartPointInIndex[1];
+
+ endPoint[0] = m_EndPointInIndex[0];
+ endPoint[1] = m_EndPointInIndex[1];
+
+ //minimum value in each direction for startRegion
+ IndexType startRegion;
+ startRegion[0] = startPoint[0] < endPoint[0] ? startPoint[0] : endPoint[0];
+ startRegion[1] = startPoint[1] < endPoint[1] ? startPoint[1] : endPoint[1];
+
+ //maximum value in each direction for size
+ typename InputImageType::SizeType size;
+ size[0] = abs( startPoint[0] - endPoint[0] );
+ size[1] = abs( startPoint[1] - endPoint[1] );
+
+
+ typename CostFunctionType::RegionType region;
+ region.SetSize( size );
+ region.SetIndex( startRegion );
+ /*---------------------------------------------*/
+
+ //inputImage->SetRequestedRegion(region);
+
+ typedef itk::CastImageFilter< InputImageType, FloatImageType > CastFilterType;
+ typename CastFilterType::Pointer castFilter = CastFilterType::New();
+ castFilter->SetInput(inputImage);
+ castFilter->Update();
+ /* extracts features from image and calculates costs */
+ if( m_ImageModified )
+ m_CostFunction->SetImage(castFilter->GetOutput());
+ m_CostFunction->SetStartIndex(startPoint);
+ m_CostFunction->SetEndIndex(endPoint);
+ m_CostFunction->SetRequestedRegion(region);
+ m_CostFunction->SetUseCostMap(m_UseDynamicCostMap);
+ /*---------------------------------------------*/
+
+
+ /* calculate shortest path between start and end point */
+ m_ShortestPathFilter->SetFullNeighborsMode(true);
+ m_ShortestPathFilter->SetInput(castFilter->GetOutput());
+ m_ShortestPathFilter->SetMakeOutputImage(false);
+
+ //m_ShortestPathFilter->SetCalcAllDistances(true);
+ m_ShortestPathFilter->SetStartIndex(startPoint);
+ m_ShortestPathFilter->SetEndIndex(endPoint);
+
+
+ m_ShortestPathFilter->Update();
+
+ /*---------------------------------------------*/
+
+
+ /* construct contour from path image */
+ //get the shortest path as vector
+ typename std::vector< ShortestPathImageFilterType::IndexType> shortestPath = m_ShortestPathFilter->GetVectorPath();
+
+ //fill the output contour with controll points from the path
+ OutputType::Pointer output = dynamic_cast ( this->MakeOutput( 0 ).GetPointer() );
+ this->SetNthOutput(0, output.GetPointer());
+
+ output->Expand(m_Timestep+1);
+
+ mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+
+ typename std::vector< ShortestPathImageFilterType::IndexType>::iterator pathIterator = shortestPath.begin();
+
+ while(pathIterator != shortestPath.end())
+ {
+ mitk::Point3D currentPoint;
+ currentPoint[0] = (*pathIterator)[0];
+ currentPoint[1] = (*pathIterator)[1];
+ currentPoint[2] = 0;
+
+
+ input->GetGeometry()->IndexToWorld(currentPoint, currentPoint);
+ output->AddVertex(currentPoint, false, m_Timestep);
+
+ pathIterator++;
+ }
+ /*---------------------------------------------*/
+}
+
+bool mitk::ImageLiveWireContourModelFilter::CreateDynamicCostMap(mitk::ContourModel* path)
+{
+ mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+ if(input)
+ {
+ AccessFixedDimensionByItk_1(input,CreateDynamicCostMapByITK, 2, path);
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+
+template
+void mitk::ImageLiveWireContourModelFilter::CreateDynamicCostMapByITK( itk::Image* inputImage, mitk::ContourModel* path )
+{
+ /*++++++++++ create dynamic cost transfer map ++++++++++*/
+
+ /* Compute the costs of the gradient magnitude dynamically.
+ * using a map of the histogram of gradient magnitude image.
+ * Use the histogram gradient map to interpolate the costs
+ * with gaussing function including next two bins right and left
+ * to current position x. With the histogram gradient costs are interpolated
+ * with a gaussing function summation of next two bins right and left
+ * to current position x.
+ */
+ std::vector< itk::Index > shortestPath;
+
+ mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+ if(path == NULL)
+ {
+ OutputType::Pointer output = this->GetOutput();
+ mitk::ContourModel::VertexIterator it = output->IteratorBegin();
+ while( it != output->IteratorEnd() )
+ {
+ itk::Index cur;
+ mitk::Point3D c = (*it)->Coordinates;
+ input->GetGeometry()->WorldToIndex(c, c);
+ cur[0] = c[0];
+ cur[1] = c[1];
+
+ shortestPath.push_back( cur);
+ it++;
+ }
+ }
+ else
+ {
+
+ mitk::ContourModel::VertexIterator it = path->IteratorBegin();
+ while( it != path->IteratorEnd() )
+ {
+ itk::Index cur;
+ mitk::Point3D c = (*it)->Coordinates;
+ input->GetGeometry()->WorldToIndex(c, c);
+ cur[0] = c[0];
+ cur[1] = c[1];
+
+ shortestPath.push_back( cur);
+ it++;
+ }
+
+ }
+
+
+ /*+++ filter image gradient magnitude +++*/
+ typedef itk::GradientMagnitudeImageFilter< itk::Image, itk::Image > GradientMagnitudeFilterType;
+ typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New();
+ gradientFilter->SetInput(inputImage);
+ gradientFilter->Update();
+ typename itk::Image::Pointer gradientMagnImage = gradientFilter->GetOutput();
+
+ //get the path
+
+
+ //iterator of path
+ typename std::vector< itk::Index >::iterator pathIterator = shortestPath.begin();
+
+ std::map< int, int > histogram;
+
+ //create histogram within path
+ while(pathIterator != shortestPath.end())
+ {
+ //count pixel values
+ //use scale factor to avoid mapping gradients between 0.0 and 1.0 to same bin
+ histogram[ static_cast( gradientMagnImage->GetPixel((*pathIterator)) * ImageLiveWireContourModelFilter::CostFunctionType::MAPSCALEFACTOR ) ] += 1;
+
+ pathIterator++;
+ }
+
+ double max = 1.0;
+
+ if( !histogram.empty() )
+ {
+
+ std::map< int, int >::iterator itMAX;
+
+ //get max of histogramm
+ int currentMaxValue = 0;
+ std::map< int, int >::iterator it = histogram.begin();
+ while( it != histogram.end())
+ {
+ if((*it).second > currentMaxValue)
+ {
+ itMAX = it;
+ currentMaxValue = (*it).second;
+ }
+ it++;
+ }
+
+
+ std::map< int, int >::key_type keyOfMax = itMAX->first;
+
+
+ /*+++++++++++++++++++++++++ compute the to max of gaussian summation ++++++++++++++++++++++++*/
+ std::map< int, int >::iterator end = histogram.end();
+ std::map< int, int >::iterator last = --(histogram.end());
+
+ std::map< int, int >::iterator left2;
+ std::map< int, int >::iterator left1;
+ std::map< int, int >::iterator right1;
+ std::map< int, int >::iterator right2;
+
+ right1 = itMAX;
+
+
+ if(right1 == end || right1 == last )
+ {
+ right2 = end;
+ }
+ else//( right1 <= last )
+ {
+ std::map< int, int >::iterator temp = right1;
+ right2 = ++right1;//rght1 + 1
+ right1 = temp;
+ }
+
+
+ if( right1 == histogram.begin() )
+ {
+ left1 = end;
+ left2 = end;
+ }
+ else if( right1 == (++(histogram.begin())) )
+ {
+ std::map< int, int >::iterator temp = right1;
+ left1 = --right1;//rght1 - 1
+ right1 = temp;
+ left2 = end;
+ }
+ else
+ {
+ std::map< int, int >::iterator temp = right1;
+ left1 = --right1;//rght1 - 1
+ left2 = --right1;//rght1 - 2
+ right1 = temp;
+ }
+
+ double partRight1, partRight2, partLeft1, partLeft2;
+ partRight1 = partRight2 = partLeft1 = partLeft2 = 0.0;
+
+
+ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ f(x) = v(bin) * e^ ( -1/2 * (|x-k(bin)| / sigma)^2 )
+
+ gaussian approximation
+
+ where
+ v(bin) is the value in the map
+ k(bin) is the key
+ */
+ if( left2 != end )
+ {
+ partLeft2 = ImageLiveWireContourModelFilter::CostFunctionType::Gaussian(keyOfMax, left2->first, left2->second);
+ }
+
+ if( left1 != end )
+ {
+ partLeft1 = ImageLiveWireContourModelFilter::CostFunctionType::Gaussian(keyOfMax, left1->first, left1->second);
+ }
+
+ if( right1 != end )
+ {
+ partRight1 = ImageLiveWireContourModelFilter::CostFunctionType::Gaussian(keyOfMax, right1->first, right1->second);
+ }
+
+ if( right2 != end )
+ {
+ partRight2 = ImageLiveWireContourModelFilter::CostFunctionType::Gaussian(keyOfMax, right2->first, right2->second);
+ }
+ /*----------------------------------------------------------------------------*/
+
+ max = (partRight1 + partRight2 + partLeft1 + partLeft2);
+
+ }
+
+ this->m_CostFunction->SetDynamicCostMap(histogram);
+ this->m_CostFunction->SetCostMapMaximum(max);
+
+}
diff --git a/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.h b/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.h
new file mode 100644
index 0000000000..74ff4fa4e8
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageLiveWireContourModelFilter.h
@@ -0,0 +1,153 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef _mitkImageLiveWireContourModelFilter_h__
+#define _mitkImageLiveWireContourModelFilter_h__
+
+#include "mitkCommon.h"
+#include "SegmentationExports.h"
+#include "mitkContourModel.h"
+#include "mitkContourModelSource.h"
+
+#include
+#include
+#include
+
+#include
+#include
+
+
+namespace mitk {
+
+ /**
+
+ \brief Calculates a LiveWire contour between two points in an image.
+
+ For defining costs between two pixels specific features are extraced from the image and tranformed into a single cost value.
+ \sa ShortestPathCostFunctionLiveWire
+
+ The filter is able to create dynamic cost tranfer map and thus use on the fly training.
+ \Note On the fly training will only be used for next update.
+ The computation uses the last calculated segment to map cost according to features in the area of the segment.
+
+ For time resolved purposes use ImageLiveWireContourModelFilter::SetTimestep( unsigned int ) to create the LiveWire contour
+ at a specific timestep.
+
+ \ingroup ContourModelFilters
+ \ingroup Process
+ */
+ class Segmentation_EXPORT ImageLiveWireContourModelFilter : public ContourModelSource
+ {
+
+ public:
+
+ mitkClassMacro(ImageLiveWireContourModelFilter, ContourModelSource);
+ itkNewMacro(Self);
+
+
+ typedef ContourModel OutputType;
+ typedef OutputType::Pointer OutputTypePointer;
+ typedef mitk::Image InputType;
+
+ typedef itk::Image< float, 2 > FloatImageType;
+ typedef itk::ShortestPathImageFilter< FloatImageType, FloatImageType > ShortestPathImageFilterType;
+ typedef itk::ShortestPathCostFunctionLiveWire< FloatImageType > CostFunctionType;
+
+
+ /** \brief start point in worldcoordinates*/
+ itkSetMacro(StartPoint, mitk::Point3D);
+ itkGetMacro(StartPoint, mitk::Point3D);
+
+ /** \brief end point in woorldcoordinates*/
+ itkSetMacro(EndPoint, mitk::Point3D);
+ itkGetMacro(EndPoint, mitk::Point3D);
+
+ /** \brief Create dynamic cost tranfer map - use on the fly training.
+ \Note On the fly training will be used for next update only.
+ The computation uses the last calculated segment to map cost according to features in the area of the segment.
+ */
+ itkSetMacro(UseDynamicCostMap, bool);
+ itkGetMacro(UseDynamicCostMap, bool);
+
+
+ virtual void SetInput( const InputType *input);
+
+ virtual void SetInput( unsigned int idx, const InputType * input);
+
+ const InputType* GetInput(void);
+
+ const InputType* GetInput(unsigned int idx);
+
+ virtual OutputType* GetOutput();
+
+
+ /** \brief Create dynamic cost tranfer map - on the fly training*/
+ bool CreateDynamicCostMap(mitk::ContourModel* path=NULL);
+
+ void SetTimestep( unsigned int timestep )
+ {
+ m_Timestep = timestep;
+ }
+
+ unsigned int GetTimestep()
+ {
+ return m_Timestep;
+ }
+
+ protected:
+ ImageLiveWireContourModelFilter();
+
+ virtual ~ImageLiveWireContourModelFilter();
+
+ void GenerateOutputInformation() {};
+
+ void GenerateData();
+
+ /** \brief start point in worldcoordinates*/
+ mitk::Point3D m_StartPoint;
+
+ /** \brief end point in woorldcoordinates*/
+ mitk::Point3D m_EndPoint;
+
+ /** \brief Start point in index*/
+ mitk::Point3D m_StartPointInIndex;
+
+ /** \brief End point in index*/
+ mitk::Point3D m_EndPointInIndex;
+
+ /** \brief The cost function to compute costs between two pixels*/
+ CostFunctionType::Pointer m_CostFunction;
+
+ /** \brief Shortest path filter according to cost function m_CostFunction*/
+ ShortestPathImageFilterType::Pointer m_ShortestPathFilter;
+
+ /** \brief Flag to use a dynmic cost map or not*/
+ bool m_UseDynamicCostMap;
+
+ bool m_ImageModified;
+
+ unsigned int m_Timestep;
+
+ template
+ void ItkProcessImage (itk::Image* inputImage);
+
+ template
+ void CreateDynamicCostMapByITK(itk::Image* inputImage, mitk::ContourModel* path=NULL);
+ };
+
+}
+
+#endif
diff --git a/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.cpp b/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.cpp
new file mode 100644
index 0000000000..5468ca98c7
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.cpp
@@ -0,0 +1,76 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include "mitkImageToContourModelFilter.h"
+
+
+
+mitk::ImageToContourModelFilter::ImageToContourModelFilter()
+{
+}
+
+
+
+mitk::ImageToContourModelFilter::~ImageToContourModelFilter()
+{
+
+}
+
+
+
+void mitk::ImageToContourModelFilter::SetInput ( const mitk::ImageToContourModelFilter::InputType* input )
+{
+ this->SetInput( 0, input );
+}
+
+
+
+void mitk::ImageToContourModelFilter::SetInput ( unsigned int idx, const mitk::ImageToContourModelFilter::InputType* input )
+{
+ if ( idx + 1 > this->GetNumberOfInputs() )
+ {
+ this->SetNumberOfRequiredInputs(idx + 1);
+ }
+ if ( input != static_cast ( this->ProcessObject::GetInput ( idx ) ) )
+ {
+ this->ProcessObject::SetNthInput ( idx, const_cast ( input ) );
+ this->Modified();
+ }
+}
+
+
+
+const mitk::ImageToContourModelFilter::InputType* mitk::ImageToContourModelFilter::GetInput( void )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(0));
+}
+
+
+
+const mitk::ImageToContourModelFilter::InputType* mitk::ImageToContourModelFilter::GetInput( unsigned int idx )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(idx));
+}
+
+
+void mitk::ImageToContourModelFilter::GenerateData()
+{
+
+}
diff --git a/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.h b/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.h
new file mode 100644
index 0000000000..40fc5d9a4f
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageToContourModelFilter.h
@@ -0,0 +1,66 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef _mitkImageToContourModelFilter_h__
+#define _mitkImageToContourModelFilter_h__
+
+#include "mitkCommon.h"
+#include "SegmentationExports.h"
+#include "mitkContourModel.h"
+#include "mitkContourModelSource.h"
+#include
+
+
+namespace mitk {
+
+ /**
+ *
+ * \brief Base class for all filters with mitk::Image as input and mitk::ContourModel
+ *
+ * \ingroup ContourModelFilters
+ * \ingroup Process
+ */
+ class Segmentation_EXPORT ImageToContourModelFilter : public ContourModelSource
+ {
+
+ public:
+
+ mitkClassMacro(ImageToContourModelFilter, ContourModelSource);
+ itkNewMacro(Self);
+
+ typedef mitk::Image InputType;
+
+
+ virtual void SetInput( const InputType *input);
+
+ virtual void SetInput( unsigned int idx, const InputType * input);
+
+ const InputType* GetInput(void);
+
+ const InputType* GetInput(unsigned int idx);
+
+ protected:
+ ImageToContourModelFilter();
+
+ virtual ~ImageToContourModelFilter();
+
+ void GenerateData();
+
+ };
+
+}
+
+#endif
diff --git a/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.cpp b/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.cpp
new file mode 100644
index 0000000000..a9b7a18f0d
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.cpp
@@ -0,0 +1,194 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include "mitkImageToLiveWireContourFilter.h"
+
+
+#include
+#include
+#include
+
+
+
+mitk::ImageToLiveWireContourFilter::ImageToLiveWireContourFilter()
+{
+ OutputType::Pointer output = dynamic_cast ( this->MakeOutput( 0 ).GetPointer() );
+ this->SetNumberOfRequiredInputs(1);
+ this->SetNumberOfOutputs( 1 );
+ this->SetNthOutput(0, output.GetPointer());
+}
+
+
+
+mitk::ImageToLiveWireContourFilter::~ImageToLiveWireContourFilter()
+{
+
+}
+
+
+
+void mitk::ImageToLiveWireContourFilter::SetInput ( const mitk::ImageToLiveWireContourFilter::InputType* input )
+{
+ this->SetInput( 0, input );
+}
+
+void mitk::ImageToLiveWireContourFilter::SetInput ( unsigned int idx, const mitk::ImageToLiveWireContourFilter::InputType* input )
+{
+ if ( idx + 1 > this->GetNumberOfInputs() )
+ {
+ this->SetNumberOfRequiredInputs(idx + 1);
+ }
+ if ( input != static_cast ( this->ProcessObject::GetInput ( idx ) ) )
+ {
+ this->ProcessObject::SetNthInput ( idx, const_cast ( input ) );
+ this->Modified();
+ }
+}
+
+
+
+const mitk::ImageToLiveWireContourFilter::InputType* mitk::ImageToLiveWireContourFilter::GetInput( void )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(0));
+}
+
+
+const mitk::ImageToLiveWireContourFilter::InputType* mitk::ImageToLiveWireContourFilter::GetInput( unsigned int idx )
+{
+ if (this->GetNumberOfInputs() < 1)
+ return NULL;
+ return static_cast(this->ProcessObject::GetInput(idx));
+}
+
+
+
+void mitk::ImageToLiveWireContourFilter::GenerateData()
+{
+ mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+
+ if(!input)
+ {
+ MITK_ERROR << "No input available.";
+ itkExceptionMacro("mitk::ImageToLiveWireContourFilter: No input available. Please set the input!");
+ return;
+ }
+
+ if( input->GetDimension() != 2 )
+ {
+ MITK_ERROR << "Filter is only working on 2D images.";
+ itkExceptionMacro("mitk::ImageToLiveWireContourFilter: Filter is only working on 2D images.. Please make sure that the input is 2D!");
+ return;
+ }
+
+
+ input->GetGeometry()->WorldToIndex(m_StartPoint, m_StartPointInIndex);
+ input->GetGeometry()->WorldToIndex(m_EndPoint, m_EndPointInIndex);
+
+
+ AccessFixedDimensionByItk(input, ItkProcessImage, 2);
+}
+
+
+
+template
+void mitk::ImageToLiveWireContourFilter::ItkProcessImage (itk::Image* inputImage)
+{
+ //typedef itk::Image< TPixel, VImageDimension > InputImageType;
+ //typedef itk::Image< float, 2 > FloatImageType;
+
+ //typedef typename itk::ShortestPathImageFilter< InputImageType, InputImageType > ShortestPathImageFilterType;
+ //typedef typename itk::ShortestPathCostFunctionLiveWire< InputImageType > CostFunctionType;
+
+ //typedef InputImageType::IndexType IndexType;
+
+
+ ///* compute the requested region for itk filters */
+
+ //typename IndexType startPoint, endPoint;
+ //
+ //startPoint[0] = m_StartPointInIndex[0];
+ //startPoint[1] = m_StartPointInIndex[1];
+
+ //endPoint[0] = m_EndPointInIndex[0];
+ //endPoint[1] = m_EndPointInIndex[1];
+
+ ////minimum value in each direction for startRegion
+ //typename IndexType startRegion;
+ //startRegion[0] = startPoint[0] < endPoint[0] ? startPoint[0] : endPoint[0];
+ //startRegion[1] = startPoint[1] < endPoint[1] ? startPoint[1] : endPoint[1];
+
+ ////maximum value in each direction for size
+ //typename InputImageType::SizeType size;
+ //size[0] = startPoint[0] > endPoint[0] ? startPoint[0] : endPoint[0];
+ //size[1] = startPoint[1] > endPoint[1] ? startPoint[1] : endPoint[1];
+
+
+ //typename InputImageType::RegionType region;
+ //region.SetSize( size );
+ //region.SetIndex( startRegion );
+ ///*---------------------------------------------*/
+
+
+ ///* extracts features from image and calculates costs */
+ //typename CostFunctionType::Pointer costFunction = CostFunctionType::New();
+ //costFunction->SetImage(inputImage);
+ //costFunction->SetStartIndex(startPoint);
+ //costFunction->SetEndIndex(endPoint);
+ //costFunction->SetRequestedRegion(region);
+ ///*---------------------------------------------*/
+
+
+ ///* calculate shortest path between start and end point */
+ //ShortestPathImageFilterType::Pointer shortestPathFilter = ShortestPathImageFilterType::New();
+ //shortestPathFilter->SetFullNeighborsMode(true);
+ //shortestPathFilter->SetInput(inputImage);
+ //shortestPathFilter->SetMakeOutputImage(true);
+ //shortestPathFilter->SetStoreVectorOrder(false);
+ ////shortestPathFilter->SetActivateTimeOut(true);
+ //shortestPathFilter->SetStartIndex(startPoint);
+ //shortestPathFilter->SetEndIndex(endPoint);
+
+ //shortestPathFilter->Update();
+
+ ///*---------------------------------------------*/
+
+
+ ///* construct contour from path image */
+ ////get the shortest path as vector
+ //std::vector< itk::Index<3> > shortestPath = shortestPathFilter->GetVectorPath();
+
+ ////fill the output contour with controll points from the path
+ //OutputType::Pointer outputContour = this->GetOutput();
+ //mitk::Image::ConstPointer input = dynamic_cast(this->GetInput());
+
+ //std::vector< itk::Index<3> >::iterator pathIterator = shortestPath.begin();
+
+ //while(pathIterator != shortestPath.end())
+ //{
+ // mitk::Point3D currentPoint;
+ // currentPoint[0] = (*pathIterator)[0];
+ // currentPoint[1] = (*pathIterator)[1];
+
+ // input->GetGeometry(0)->IndexToWorld(currentPoint, currentPoint);
+ // outputContour->AddVertex(currentPoint);
+ //
+ // pathIterator++;
+ //}
+ /*---------------------------------------------*/
+
+}
diff --git a/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.h b/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.h
new file mode 100644
index 0000000000..eebf622d6b
--- /dev/null
+++ b/Modules/Segmentation/Algorithms/mitkImageToLiveWireContourFilter.h
@@ -0,0 +1,92 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef _mitkImageToLiveWireContourFilter_h__
+#define _mitkImageToLiveWireContourFilter_h__
+
+#include "mitkCommon.h"
+#include "SegmentationExports.h"
+#include "mitkContourModel.h"
+#include "mitkContourModelSource.h"
+
+#include
+#include
+#include
+
+
+
+
+namespace mitk {
+
+ /**
+ *
+ * \brief
+ *
+ * \ingroup ContourModelFilters
+ * \ingroup Process
+ */
+ class Segmentation_EXPORT ImageToLiveWireContourFilter : public ContourModelSource
+ {
+
+ public:
+
+ mitkClassMacro(ImageToLiveWireContourFilter, ContourModelSource);
+ itkNewMacro(Self);
+
+ typedef ContourModel OutputType;
+ typedef OutputType::Pointer OutputTypePointer;
+ typedef mitk::Image InputType;
+
+ itkSetMacro(StartPoint, mitk::Point3D);
+ itkGetMacro(StartPoint, mitk::Point3D);
+
+ itkSetMacro(EndPoint, mitk::Point3D);
+ itkGetMacro(EndPoint, mitk::Point3D);
+
+
+ virtual void SetInput( const InputType *input);
+
+ virtual void SetInput( unsigned int idx, const InputType * input);
+
+ const InputType* GetInput(void);
+
+ const InputType* GetInput(unsigned int idx);
+
+ protected:
+ ImageToLiveWireContourFilter();
+
+ virtual ~ImageToLiveWireContourFilter();
+
+ void GenerateData();
+
+ void GenerateOutputInformation() {};
+
+ mitk::Point3D m_StartPoint;
+ mitk::Point3D m_EndPoint;
+
+ mitk::Point3D m_StartPointInIndex;
+ mitk::Point3D m_EndPointInIndex;
+
+ private:
+
+ template
+ void ItkProcessImage (itk::Image* inputImage);
+
+
+ };
+
+}
+#endif
diff --git a/Modules/Segmentation/CMakeLists.txt b/Modules/Segmentation/CMakeLists.txt
index 133b39a60c..f156909325 100644
--- a/Modules/Segmentation/CMakeLists.txt
+++ b/Modules/Segmentation/CMakeLists.txt
@@ -1,10 +1,10 @@
#configure_file(${PROJECT_SOURCE_DIR}/CMake/ToolExtensionITKFactory.cpp.in $#{PROJECT_BINARY_DIR}/ToolExtensionITKFactory.cpp.in COPYONLY)
#configure_file(${PROJECT_SOURCE_DIR}/CMake/ToolExtensionITKFactoryLoader.cpp.in $#{PROJECT_BINARY_DIR}/ToolExtensionITKFactoryLoader.cpp.in COPYONLY)
#configure_file(${PROJECT_SOURCE_DIR}/CMake/ToolGUIExtensionITKFactory.cpp.in $#{PROJECT_BINARY_DIR}/ToolGUIExtensionITKFactory.cpp.in COPYONLY)
MITK_CREATE_MODULE( Segmentation
INCLUDE_DIRS Algorithms Controllers DataManagement Interactions IO Rendering
- DEPENDS Mitk ipSegmentation mitkIpFunc MitkExt
+ DEPENDS Mitk ipSegmentation mitkIpFunc MitkExt MitkGraphAlgorithms
)
add_subdirectory(Testing)
\ No newline at end of file
diff --git a/Modules/Segmentation/DataManagement/mitkContourElement.cpp b/Modules/Segmentation/DataManagement/mitkContourElement.cpp
index 02357705ba..f617e988ba 100644
--- a/Modules/Segmentation/DataManagement/mitkContourElement.cpp
+++ b/Modules/Segmentation/DataManagement/mitkContourElement.cpp
@@ -1,351 +1,356 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include
mitk::ContourElement::ContourElement()
{
this->m_Vertices = new VertexListType();
this->m_IsClosed = false;
}
mitk::ContourElement::ContourElement(const mitk::ContourElement &other) :
m_Vertices(other.m_Vertices), m_IsClosed(other.m_IsClosed)
{
}
mitk::ContourElement::~ContourElement()
{
delete this->m_Vertices;
}
void mitk::ContourElement::AddVertex(mitk::Point3D &vertex, bool isControlPoint)
{
this->m_Vertices->push_back(new VertexType(vertex, isControlPoint));
}
void mitk::ContourElement::AddVertex(VertexType &vertex)
{
this->m_Vertices->push_back(&vertex);
}
void mitk::ContourElement::AddVertexAtFront(mitk::Point3D &vertex, bool isControlPoint)
{
this->m_Vertices->push_front(new VertexType(vertex, isControlPoint));
}
void mitk::ContourElement::AddVertexAtFront(VertexType &vertex)
{
this->m_Vertices->push_front(&vertex);
}
void mitk::ContourElement::InsertVertexAtIndex(mitk::Point3D &vertex, bool isControlPoint, int index)
{
if(index > 0 && this->GetSize() > index)
{
VertexIterator _where = this->m_Vertices->begin();
_where += index;
this->m_Vertices->insert(_where, new VertexType(vertex, isControlPoint));
}
}
mitk::ContourElement::VertexType* mitk::ContourElement::GetVertexAt(int index)
{
return this->m_Vertices->at(index);
}
mitk::ContourElement::VertexType* mitk::ContourElement::GetVertexAt(const mitk::Point3D &point, float eps)
{
/* current version iterates over the whole deque - should some kind of an octree with spatial query*/
if(eps > 0)
{
+ if(this->m_Vertices->size() == 1)
+ {
+ return this->m_Vertices->at(0);
+ }
+
if(true ) //currently no method with better performance is available
{
return BruteForceGetVertexAt(point, eps);
}
else
{
return OptimizedGetVertexAt(point, eps);
}
}//if eps < 0
return NULL;
}
mitk::ContourElement::VertexType* mitk::ContourElement::BruteForceGetVertexAt(const mitk::Point3D &point, float eps)
{
if(eps > 0)
{
std::deque< std::pair > nearestlist;
ConstVertexIterator it = this->m_Vertices->begin();
ConstVertexIterator end = this->m_Vertices->end();
while(it != end)
{
mitk::Point3D currentPoint = (*it)->Coordinates;
double distance = currentPoint.EuclideanDistanceTo(point);
if(distance < eps)
{
//if list is emtpy, add point to list
if(nearestlist.size() < 1)
{
nearestlist.push_front(std::pair( (*it)->Coordinates.EuclideanDistanceTo(point), (*it) ));
}
//found an approximate point - check if current is closer then first in nearestlist
else if( distance < nearestlist.front().first )
{
//found even closer vertex
nearestlist.push_front(std::pair( (*it)->Coordinates.EuclideanDistanceTo(point), (*it) ));
}
}//if distance > eps
it++;
}//while
if(nearestlist.size() > 0)
{
/*++++++++++++++++++++ return the nearest active point if one was found++++++++++++++++++*/
std::deque< std::pair >::iterator it = nearestlist.begin();
std::deque< std::pair >::iterator end = nearestlist.end();
while(it != end)
{
if( (*it).second->IsControlPoint )
{
return (*it).second;
}
it++;
}
/*---------------------------------------------------------------------------------------*/
//return closest point
return nearestlist.front().second;
}
}
return NULL;
}
mitk::ContourElement::VertexType* mitk::ContourElement::OptimizedGetVertexAt(const mitk::Point3D &point, float eps)
{
if( (eps > 0) && (this->m_Vertices->size()>0) )
{
int k = 1;
int dim = 3;
int nPoints = this->m_Vertices->size();
ANNpointArray pointsArray;
ANNpoint queryPoint;
ANNidxArray indexArray;
ANNdistArray distanceArray;
ANNkd_tree* kdTree;
queryPoint = annAllocPt(dim);
pointsArray = annAllocPts(nPoints, dim);
indexArray = new ANNidx[k];
distanceArray = new ANNdist[k];
int i = 0;
//fill points array with our control points
for(VertexIterator it = this->m_Vertices->begin(); it != this->m_Vertices->end(); it++, i++)
{
mitk::Point3D cur = (*it)->Coordinates;
pointsArray[i][0]= cur[0];
pointsArray[i][1]= cur[1];
pointsArray[i][2]= cur[2];
}
//create the kd tree
kdTree = new ANNkd_tree(pointsArray,nPoints, dim);
//fill mitk::Point3D into ANN query point
queryPoint[0] = point[0];
queryPoint[1] = point[1];
queryPoint[2] = point[2];
//k nearest neighbour search
kdTree->annkSearch(queryPoint, k, indexArray, distanceArray, eps);
VertexType* ret = NULL;
try
{
ret = this->m_Vertices->at(indexArray[0]);
}
catch(std::out_of_range ex)
{
//ret stays NULL
return ret;
}
//clean up ANN
delete [] indexArray;
delete [] distanceArray;
delete kdTree;
annClose();
return ret;
}
return NULL;
}
mitk::ContourElement::VertexListType* mitk::ContourElement::GetVertexList()
{
return this->m_Vertices;
}
bool mitk::ContourElement::IsClosed()
{
return this->m_IsClosed;
}
void mitk::ContourElement::Close()
{
this->m_IsClosed = true;
}
void mitk::ContourElement::Open()
{
this->m_IsClosed = false;
}
void mitk::ContourElement::SetIsClosed( bool isClosed)
{
isClosed ? this->Close() : this->Open();
}
void mitk::ContourElement::Concatenate(mitk::ContourElement* other)
{
if( other->GetSize() > 0)
{
ConstVertexIterator it = other->m_Vertices->begin();
ConstVertexIterator end = other->m_Vertices->end();
//add all vertices of other after last vertex
while(it != end)
{
this->m_Vertices->push_back(*it);
it++;
}
}
}
bool mitk::ContourElement::RemoveVertex(mitk::ContourElement::VertexType* vertex)
{
VertexIterator it = this->m_Vertices->begin();
VertexIterator end = this->m_Vertices->end();
//search for vertex and remove it if exists
while(it != end)
{
if((*it) == vertex)
{
this->m_Vertices->erase(it);
return true;
}
it++;
}
return false;
}
bool mitk::ContourElement::RemoveVertexAt(int index)
{
if( index >= 0 && index < this->m_Vertices->size() )
{
this->m_Vertices->erase(this->m_Vertices->begin()+index);
return true;
}
else
{
return false;
}
}
bool mitk::ContourElement::RemoveVertexAt(mitk::Point3D &point, float eps)
{
/* current version iterates over the whole deque - should be some kind of an octree with spatial query*/
if(eps > 0){
VertexIterator it = this->m_Vertices->begin();
VertexIterator end = this->m_Vertices->end();
while(it != end)
{
mitk::Point3D currentPoint = (*it)->Coordinates;
if(currentPoint.EuclideanDistanceTo(point) < eps)
{
//approximate point found
//now erase it
this->m_Vertices->erase(it);
return true;
}
it++;
}
}
return false;
}
void mitk::ContourElement::Clear()
{
this->m_Vertices->clear();
}
diff --git a/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.cpp b/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.cpp
new file mode 100644
index 0000000000..71fdceeff9
--- /dev/null
+++ b/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.cpp
@@ -0,0 +1,349 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include "mitkContourModelLiveWireInteractor.h"
+
+#include "mitkToolManager.h"
+
+#include "mitkBaseRenderer.h"
+#include "mitkRenderingManager.h"
+
+#include
+#include
+
+#include
+
+
+mitk::ContourModelLiveWireInteractor::ContourModelLiveWireInteractor(DataNode* dataNode)
+:ContourModelInteractor(dataNode)
+{
+ m_LiveWireFilter = mitk::ImageLiveWireContourModelFilter::New();
+
+ mitk::ContourModel::Pointer m_ContourLeft = mitk::ContourModel::New();
+ mitk::ContourModel::Pointer m_ContourRight = mitk::ContourModel::New();
+ m_NextActiveVertexDown.Fill(0);
+ m_NextActiveVertexUp.Fill(0);
+}
+
+
+mitk::ContourModelLiveWireInteractor::~ContourModelLiveWireInteractor()
+{
+}
+
+
+
+bool mitk::ContourModelLiveWireInteractor::OnCheckPointClick( Action* action, const StateEvent* stateEvent)
+{
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) {
+ this->HandleEvent( new mitk::StateEvent(EIDNO, stateEvent->GetEvent()) );
+ return false;
+ }
+
+
+
+ mitk::StateEvent* newStateEvent = NULL;
+
+ int timestep = stateEvent->GetEvent()->GetSender()->GetTimeStep();
+
+ mitk::ContourModel *contour = dynamic_cast(
+ m_DataNode->GetData() );
+
+
+ contour->Deselect();
+
+ /*
+ * Check distance to any vertex.
+ * Transition YES if click close to a vertex
+ */
+ mitk::Point3D click = positionEvent->GetWorldPosition();
+
+
+ if (contour->SelectVertexAt(click, 1.5, timestep) )
+ {
+ contour->SetSelectedVertexAsControlPoint();
+
+ assert( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+ newStateEvent = new mitk::StateEvent(EIDYES, stateEvent->GetEvent());
+ m_lastMousePosition = click;
+
+
+ mitk::ContourModel *contour = dynamic_cast( m_DataNode->GetData() );
+
+ //
+ m_ContourLeft = mitk::ContourModel::New();
+
+ //get coordinates of next active vertex downwards from selected vertex
+ int downIndex = SplitContourFromSelectedVertex( contour, m_ContourLeft, false, timestep);
+
+ mitk::ContourModel::VertexIterator itDown = contour->IteratorBegin() + downIndex;
+ m_NextActiveVertexDown = (*itDown)->Coordinates;
+
+
+ //
+ m_ContourRight = mitk::ContourModel::New();
+
+ //get coordinates of next active vertex upwards from selected vertex
+ int upIndex = SplitContourFromSelectedVertex( contour, m_ContourRight, true, timestep);
+
+ mitk::ContourModel::VertexIterator itUp = contour->IteratorBegin() + upIndex;
+ m_NextActiveVertexUp = (*itUp)->Coordinates;
+
+ }
+ else
+ {
+ newStateEvent = new mitk::StateEvent(EIDNO, stateEvent->GetEvent());
+ }
+
+ this->HandleEvent( newStateEvent );
+
+ return true;
+}
+
+
+
+bool mitk::ContourModelLiveWireInteractor::OnDeletePoint( Action* action, const StateEvent* stateEvent)
+{
+
+ int timestep = stateEvent->GetEvent()->GetSender()->GetTimeStep();
+
+ mitk::ContourModel *contour = dynamic_cast( m_DataNode->GetData() );
+
+ if(contour->GetSelectedVertex())
+ {
+
+ mitk::ContourModel::Pointer newContour = mitk::ContourModel::New();
+ newContour->Expand(contour->GetTimeSteps());
+ newContour->Concatenate( m_ContourLeft, timestep );
+
+
+ //recompute contour between neighbored two active control points
+ this->m_LiveWireFilter->SetStartPoint( m_NextActiveVertexDown );
+ this->m_LiveWireFilter->SetEndPoint( m_NextActiveVertexUp );
+ this->m_LiveWireFilter->Update();
+
+ mitk::ContourModel::Pointer liveWireContour = mitk::ContourModel::New();
+ liveWireContour = this->m_LiveWireFilter->GetOutput();
+
+
+ //insert new liveWire computed points
+ newContour->Concatenate( liveWireContour, timestep );
+
+ newContour->Concatenate( m_ContourRight, timestep );
+
+ newContour->SetIsClosed(contour->IsClosed(timestep), timestep);
+
+ m_DataNode->SetData(newContour);
+
+
+ assert( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+
+ }
+ return true;
+}
+
+
+
+
+bool mitk::ContourModelLiveWireInteractor::OnMovePoint( Action* action, const StateEvent* stateEvent)
+{
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) return false;
+
+ int timestep = stateEvent->GetEvent()->GetSender()->GetTimeStep();
+
+ mitk::ContourModel *contour = dynamic_cast( m_DataNode->GetData() );
+
+ mitk::ContourModel::Pointer newContour = mitk::ContourModel::New();
+ newContour->Expand(contour->GetTimeSteps());
+
+ /*++++++++++++++ concatenate LEFT ++++++++++++++++++++++++++++++*/
+ newContour->Concatenate( m_ContourLeft, timestep );
+
+
+ mitk::Point3D currentPosition = positionEvent->GetWorldPosition();
+
+
+ /*+++++++++++++++ start computation ++++++++++++++++++++*/
+ //recompute contour between previous active vertex and selected vertex
+ this->m_LiveWireFilter->SetStartPoint( m_NextActiveVertexDown );
+ this->m_LiveWireFilter->SetEndPoint( currentPosition );
+ this->m_LiveWireFilter->Update();
+
+ mitk::ContourModel::Pointer liveWireContour = mitk::ContourModel::New();
+ liveWireContour = this->m_LiveWireFilter->GetOutput();
+
+ //remove point at current position because it is included in next livewire segment too.
+ liveWireContour->RemoveVertexAt( liveWireContour->GetNumberOfVertices(timestep) - 1, timestep);
+
+ /*++++++++++++++ concatenate LIVEWIRE ++++++++++++++++++++++++++++++*/
+ //insert new liveWire computed points
+ newContour->Concatenate( liveWireContour, timestep );
+
+
+ //recompute contour between selected vertex and next active vertex
+ this->m_LiveWireFilter->SetStartPoint( currentPosition );
+ this->m_LiveWireFilter->SetEndPoint( m_NextActiveVertexUp );
+ this->m_LiveWireFilter->Update();
+
+ liveWireContour = this->m_LiveWireFilter->GetOutput();
+
+ //make point at mouse position active again, so it is drawn
+ const_cast( liveWireContour->GetVertexAt(0, timestep) )->IsControlPoint = true;
+
+ /*++++++++++++++ concatenate RIGHT ++++++++++++++++++++++++++++++*/
+ //insert new liveWire computed points
+ newContour->Concatenate( liveWireContour, timestep );
+ /*------------------------------------------------*/
+
+
+ newContour->Concatenate( m_ContourRight, timestep );
+
+ newContour->SetIsClosed(contour->IsClosed(timestep), timestep);
+ newContour->Deselect();//just to make sure
+ m_DataNode->SetData(newContour);
+
+ this->m_lastMousePosition = positionEvent->GetWorldPosition();
+
+ assert( positionEvent->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() );
+
+ return true;
+}
+
+
+
+int mitk::ContourModelLiveWireInteractor::SplitContourFromSelectedVertex(mitk::ContourModel* sourceContour, mitk::ContourModel* destinationContour, bool fromSelectedUpwards, int timestep)
+{
+
+ mitk::ContourModel::VertexIterator end =sourceContour->IteratorEnd();
+ mitk::ContourModel::VertexIterator begin =sourceContour->IteratorBegin();
+
+ //search next active control point to left and rigth and set as start and end point for filter
+ mitk::ContourModel::VertexIterator itSelected = sourceContour->IteratorBegin();
+
+
+ //move iterator to position
+ while((*itSelected) != sourceContour->GetSelectedVertex())
+ {
+ itSelected++;
+ }
+
+ //CASE search upwards for next control point
+ if(fromSelectedUpwards)
+ {
+ mitk::ContourModel::VertexIterator itUp = itSelected;
+
+ if(itUp != end)
+ {
+ itUp++;//step once up otherwise the the loop breaks immediately
+ }
+
+ while( itUp != end && !((*itUp)->IsControlPoint)){ itUp++; }
+
+ mitk::ContourModel::VertexIterator it = itUp;
+
+
+ if(itSelected != begin)
+ {
+ //copy the rest of the original contour
+ while(it != end)
+ {
+ destinationContour->AddVertex( (*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ it++;
+ }
+ }
+ //else do not copy the contour
+
+
+ //return the offset of iterator at one before next-vertex-upwards
+ if(itUp != begin)
+ {
+ return std::distance( begin, itUp) - 1;
+ }
+ else
+ {
+ return std::distance( begin, itUp);
+ }
+
+ }
+ else //CASE search downwards for next control point
+ {
+ mitk::ContourModel::VertexIterator itDown = itSelected;
+ mitk::ContourModel::VertexIterator it = sourceContour->IteratorBegin();
+
+ if( itSelected != begin )
+ {
+ if(itDown != begin)
+ {
+ itDown--;//step once down otherwise the the loop breaks immediately
+ }
+
+ while( itDown != begin && !((*itDown)->IsControlPoint)){ itDown--; }
+
+ if(it != end)//if not empty
+ {
+ //always add the first vertex
+ destinationContour->AddVertex( (*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ it++;
+ }
+ //copy from begin to itDown
+ while(it <= itDown)
+ {
+ destinationContour->AddVertex( (*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ it++;
+ }
+ }
+ else
+ {
+ //if selected vertex is the first element search from end of contour downwards
+ itDown = end;
+ itDown--;
+ while(!((*itDown)->IsControlPoint)){itDown--;}
+
+ //move one forward as we don't want the first control point
+ it++;
+ //move iterator to second control point
+ while( (it!=end) && !((*it)->IsControlPoint) ){it++;}
+ //copy from begin to itDown
+ while(it <= itDown)
+ {
+ //copy the contour from second control point to itDown
+ destinationContour->AddVertex( (*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ it++;
+ }
+ }
+
+
+ //add vertex at itDown - it's not considered during while loop
+ if( it != begin && it != end)
+ {
+ //destinationContour->AddVertex( (*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ }
+
+ //return the offset of iterator at one after next-vertex-downwards
+ if( itDown != end)
+ {
+ return std::distance( begin, itDown);// + 1;//index of next vertex
+ }
+ else
+ {
+ return std::distance( begin, itDown) - 1;
+ }
+ }
+}
diff --git a/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.h b/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.h
new file mode 100644
index 0000000000..ae44899f95
--- /dev/null
+++ b/Modules/Segmentation/Interactions/mitkContourModelLiveWireInteractor.h
@@ -0,0 +1,89 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef mitkContourModelLiveWireInteractor_h_Included
+#define mitkContourModelLiveWireInteractor_h_Included
+
+#include "mitkCommon.h"
+#include "SegmentationExports.h"
+#include "mitkContourModelInteractor.h"
+#include
+#include
+
+#include
+
+namespace mitk
+{
+
+
+ /**
+ \brief
+
+ \sa Interactor
+ \sa ContourModelInteractor
+
+ \ingroup Interaction
+
+
+ \warning Make sure the working image is properly set, otherwise the algorithm for computing livewire contour segments will not work!
+
+ */
+ class Segmentation_EXPORT ContourModelLiveWireInteractor : public ContourModelInteractor
+ {
+ public:
+
+ mitkClassMacro(ContourModelLiveWireInteractor, ContourModelInteractor);
+ mitkNewMacro1Param(Self, DataNode*);
+
+
+ virtual void SetWorkingImage (mitk::Image* _arg)
+ {
+ if (this->m_WorkingImage != _arg)
+ {
+ this->m_WorkingImage = _arg;
+ this->m_LiveWireFilter->SetInput(this->m_WorkingImage);
+ this->Modified();
+ }
+ }
+
+ protected:
+
+ ContourModelLiveWireInteractor(DataNode* dataNode);
+ virtual ~ContourModelLiveWireInteractor();
+
+
+ virtual bool OnDeletePoint(Action*, const StateEvent*);
+ virtual bool OnMovePoint(Action*, const StateEvent*);
+ virtual bool OnCheckPointClick( Action* action, const StateEvent* stateEvent);
+
+ int SplitContourFromSelectedVertex(mitk::ContourModel* sourceContour,
+ mitk::ContourModel* destinationContour,
+ bool fromSelectedUpwards,
+ int timestep);
+
+ mitk::ImageLiveWireContourModelFilter::Pointer m_LiveWireFilter;
+ mitk::Image::Pointer m_WorkingImage;
+
+ mitk::Point3D m_NextActiveVertexDown;
+ mitk::Point3D m_NextActiveVertexUp;
+ mitk::ContourModel::Pointer m_ContourLeft;
+ mitk::ContourModel::Pointer m_ContourRight;
+
+ };
+
+} // namespace
+
+#endif
diff --git a/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp
new file mode 100644
index 0000000000..20082f9f8e
--- /dev/null
+++ b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.cpp
@@ -0,0 +1,653 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include "mitkLiveWireTool2D.h"
+
+#include "mitkToolManager.h"
+#include "mitkBaseRenderer.h"
+#include "mitkRenderingManager.h"
+
+#include "mitkLiveWireTool2D.xpm"
+
+#include
+#include
+#include
+
+#include "mitkContourUtils.h"
+#include "mitkContour.h"
+
+#include
+
+
+namespace mitk {
+ MITK_TOOL_MACRO(Segmentation_EXPORT, LiveWireTool2D, "LiveWire tool");
+}
+
+
+
+mitk::LiveWireTool2D::LiveWireTool2D()
+:SegTool2D("LiveWireTool")
+{
+ m_Contour = mitk::ContourModel::New();
+ m_ContourModelNode = mitk::DataNode::New();
+ m_ContourModelNode->SetData( m_Contour );
+ m_ContourModelNode->SetProperty("name", StringProperty::New("working contour node"));
+ m_ContourModelNode->SetProperty("visible", BoolProperty::New(true));
+ m_ContourModelNode->AddProperty( "contour.color", ColorProperty::New(0.9, 1.0, 0.1), NULL, true );
+ m_ContourModelNode->AddProperty( "selectedcolor", ColorProperty::New(1.0, 0.0, 0.1), NULL, true );
+
+
+ m_LiveWireContour = mitk::ContourModel::New();
+ m_LiveWireContourNode = mitk::DataNode::New();
+ //m_LiveWireContourNode->SetData( m_LiveWireContour );
+ m_LiveWireContourNode->SetProperty("name", StringProperty::New("active livewire node"));
+ m_LiveWireContourNode->SetProperty("visible", BoolProperty::New(true));
+ m_LiveWireContourNode->AddProperty( "contour.color", ColorProperty::New(0.1, 1.0, 0.1), NULL, true );
+ m_LiveWireContourNode->AddProperty( "selectedcolor", ColorProperty::New(0.5, 0.5, 0.1), NULL, true );
+
+
+ m_LiveWireFilter = mitk::ImageLiveWireContourModelFilter::New();
+
+
+ // great magic numbers
+ CONNECT_ACTION( AcINITNEWOBJECT, OnInitLiveWire );
+ CONNECT_ACTION( AcADDPOINT, OnAddPoint );
+ CONNECT_ACTION( AcMOVE, OnMouseMoveNoDynamicCosts );
+ CONNECT_ACTION( AcCHECKPOINT, OnCheckPoint );
+ CONNECT_ACTION( AcFINISH, OnFinish );
+ CONNECT_ACTION( AcDELETEPOINT, OnLastSegmentDelete );
+ CONNECT_ACTION( AcADDLINE, OnMouseMoved );
+}
+
+
+
+
+mitk::LiveWireTool2D::~LiveWireTool2D()
+{
+ m_Contours.clear();
+}
+
+
+
+float mitk::LiveWireTool2D::CanHandleEvent( StateEvent const *stateEvent) const
+{
+ mitk::PositionEvent const *positionEvent =
+ dynamic_cast (stateEvent->GetEvent());
+
+ //Key event handling:
+ if (positionEvent == NULL)
+ {
+ //check for delete and escape event
+ if(stateEvent->GetId() == 12 || stateEvent->GetId() == 14)
+ {
+ return 1.0;
+ }
+ //check, if the current state has a transition waiting for that key event.
+ else if (this->GetCurrentState()->GetTransition(stateEvent->GetId())!=NULL)
+ {
+ return 0.5;
+ }
+ else
+ {
+ return 0.0;
+ }
+ }
+ else
+ {
+ if ( positionEvent->GetSender()->GetMapperID() != BaseRenderer::Standard2D )
+ return 0.0; // we don't want anything but 2D
+
+ return 1.0;
+ }
+
+}
+
+
+
+
+const char** mitk::LiveWireTool2D::GetXPM() const
+{
+ return mitkLiveWireTool2D_xpm;
+}
+
+
+
+
+const char* mitk::LiveWireTool2D::GetName() const
+{
+ return "LiveWire";
+}
+
+
+
+
+void mitk::LiveWireTool2D::Activated()
+{
+ Superclass::Activated();
+}
+
+
+
+
+void mitk::LiveWireTool2D::Deactivated()
+{
+ this->FinishTool();
+
+ DataNode* workingNode( m_ToolManager->GetWorkingData(0) );
+ if ( !workingNode ) return;
+
+ Image* workingImage = dynamic_cast(workingNode->GetData());
+ if ( !workingImage ) return;
+
+ ContourUtils::Pointer contourUtils = mitk::ContourUtils::New();
+
+ /*+++++++++++++++++++++++ for all contours in list (currently created by tool) ++++++++++++++++++++++++++++++++++++*/
+ std::vector< std::pair >::iterator it = m_Contours.begin();
+ while(it != m_Contours.end() )
+ {
+
+ //++++++++++if node contains data
+ if( it->first->GetData() )
+ {
+
+ //+++++++++++++++if this is a contourModel
+ mitk::ContourModel* contourModel = dynamic_cast(it->first->GetData());
+ if( contourModel )
+ {
+
+ //++++++++++++++++++++++ for each timestep of this contourModel
+ for( int currentTimestep = 0; currentTimestep < contourModel->GetTimeSlicedGeometry()->GetTimeSteps(); currentTimestep++)
+ {
+
+ //get the segmentation image slice at current timestep
+ mitk::Image::Pointer workingSlice = this->GetAffectedImageSliceAs2DImage(it->second, workingImage, currentTimestep);
+
+ /*++++++++++++++++++++++ transfer to plain old contour to use contour util functionality +++++++++++++++++++++++*/
+ mitk::Contour::Pointer plainOldContour = mitk::Contour::New();
+ mitk::ContourModel::VertexIterator iter = contourModel->IteratorBegin(currentTimestep);
+ while(iter != contourModel->IteratorEnd(currentTimestep) )
+ {
+ plainOldContour->AddVertex( (*iter)->Coordinates );
+ iter++;
+ }
+ /*-------------------------------------------------------------------------------*/
+
+
+ mitk::Contour::Pointer projectedContour = contourUtils->ProjectContourTo2DSlice(workingSlice, plainOldContour, true, false);
+ contourUtils->FillContourInSlice(projectedContour, workingSlice, 1.0);
+
+ //write back to image volume
+ this->WriteBackSegmentationResult(it->second, workingSlice, currentTimestep);
+ }
+
+ //remove contour node from datastorage
+ m_ToolManager->GetDataStorage()->Remove( it->first );
+ }
+ }
+
+ ++it;
+ }
+
+ m_Contours.clear();
+
+ Superclass::Deactivated();
+}
+
+
+
+
+bool mitk::LiveWireTool2D::OnInitLiveWire (Action* action, const StateEvent* stateEvent)
+{
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) return false;
+
+ m_LastEventSender = positionEvent->GetSender();
+ m_LastEventSlice = m_LastEventSender->GetSlice();
+
+ if ( Superclass::CanHandleEvent(stateEvent) < 1.0 ) return false;
+
+
+ int timestep = positionEvent->GetSender()->GetTimeStep();
+
+ m_Contour->Expand(timestep+1);
+
+ m_ToolManager->GetDataStorage()->Add( m_ContourModelNode );
+
+ m_ToolManager->GetDataStorage()->Add( m_LiveWireContourNode );
+
+ //set current slice as input for ImageToLiveWireContourFilter
+ m_WorkingSlice = this->GetAffectedReferenceSlice(positionEvent);
+ m_LiveWireFilter->SetInput(m_WorkingSlice);
+
+ //map click to pixel coordinates
+ mitk::Point3D click = const_cast(positionEvent->GetWorldPosition());
+ itk::Index<3> idx;
+ m_WorkingSlice->GetGeometry()->WorldToIndex(click, idx);
+
+ /*+++++++++++++++++++++++ get the pixel the gradient in region of 5x5 ++++++++++++++++++++++++++*/
+ itk::Index<3> indexWithHighestGradient;
+ AccessFixedDimensionByItk_2(m_WorkingSlice, FindHighestGradientMagnitudeByITK, 2, idx, indexWithHighestGradient);
+ /*----------------------------------------------------------------------------------------------------------------*/
+
+ //itk::Index to mitk::Point3D
+ click[0] = indexWithHighestGradient[0];
+ click[1] = indexWithHighestGradient[1];
+ click[2] = indexWithHighestGradient[2];
+ m_WorkingSlice->GetGeometry()->IndexToWorld(click, click);
+
+ //set initial start point
+ m_Contour->AddVertex( click, true, timestep );
+ m_LiveWireFilter->SetStartPoint(click);
+
+ m_CreateAndUseDynamicCosts = true;
+
+ //render
+ assert( positionEvent->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() );
+
+ return true;
+}
+
+
+
+
+bool mitk::LiveWireTool2D::OnAddPoint (Action* action, const StateEvent* stateEvent)
+{
+ //complete LiveWire interaction for last segment
+ //add current LiveWire contour to the finished contour and reset
+ //to start new segment and computation
+
+ /* check if event can be handled */
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) return false;
+
+ if ( Superclass::CanHandleEvent(stateEvent) < 1.0 ) return false;
+ /* END check if event can be handled */
+
+
+ int timestep = positionEvent->GetSender()->GetTimeStep();
+
+ //remove duplicate first vertex, it's already contained in m_Contour
+ m_LiveWireContour->RemoveVertexAt(0, timestep);
+
+
+ /* TODO fix this hack*/
+ //set last to active added point
+ if( m_LiveWireContour->GetNumberOfVertices(timestep) > 0)
+ {
+ const_cast( m_LiveWireContour->GetVertexAt(m_LiveWireContour->GetNumberOfVertices(timestep)-1, timestep) )->IsControlPoint = true;
+ }
+
+ //merge contours
+ m_Contour->Concatenate(m_LiveWireContour, timestep);
+
+
+ //clear the livewire contour and reset the corresponding datanode
+ m_LiveWireContour->Clear(timestep);
+
+ //set new start point
+ m_LiveWireFilter->SetStartPoint(const_cast(positionEvent->GetWorldPosition()));
+
+ if( m_CreateAndUseDynamicCosts )
+ {
+ //use dynamic cost map for next update
+ m_LiveWireFilter->CreateDynamicCostMap(m_Contour);
+ m_LiveWireFilter->SetUseDynamicCostMap(true);
+ //m_CreateAndUseDynamicCosts = false;
+ }
+
+ //render
+ assert( positionEvent->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() );
+
+ return true;
+}
+
+
+
+
+bool mitk::LiveWireTool2D::OnMouseMoved( Action* action, const StateEvent* stateEvent)
+{
+ //compute LiveWire segment from last control point to current mouse position
+
+ /* check if event can be handled */
+ if ( Superclass::CanHandleEvent(stateEvent) < 1.0 ) return false;
+
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) return false;
+ /* END check if event can be handled */
+
+
+ /* actual LiveWire computation */
+ int timestep = positionEvent->GetSender()->GetTimeStep();
+
+ m_LiveWireFilter->SetEndPoint(const_cast(positionEvent->GetWorldPosition()));
+
+ m_LiveWireFilter->SetTimestep(timestep);
+ m_LiveWireFilter->Update();
+
+
+ //ContourModel::VertexType* currentVertex = const_cast(m_LiveWireContour->GetVertexAt(0));
+
+ this->m_LiveWireContour = this->m_LiveWireFilter->GetOutput();
+ this->m_LiveWireContourNode->SetData(this->m_LiveWireFilter->GetOutput());
+
+ /* END actual LiveWire computation */
+
+ //render
+ assert( positionEvent->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() );
+
+ return true;
+}
+
+
+
+bool mitk::LiveWireTool2D::OnMouseMoveNoDynamicCosts(Action* action, const StateEvent* stateEvent)
+{
+ //do not use dynamic cost map
+ m_LiveWireFilter->SetUseDynamicCostMap(false);
+ OnMouseMoved(action, stateEvent);
+ m_LiveWireFilter->SetUseDynamicCostMap(true);
+ return true;
+}
+
+
+
+
+
+bool mitk::LiveWireTool2D::OnCheckPoint( Action* action, const StateEvent* stateEvent)
+{
+ //check double click on first control point to finish the LiveWire tool
+ //
+ //Check distance to first point.
+ //Transition YES if click close to first control point
+ //
+
+ mitk::StateEvent* newStateEvent = NULL;
+
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent)
+ {
+ //stay in current state
+ newStateEvent = new mitk::StateEvent(EIDNO, stateEvent->GetEvent());
+ }
+ else
+ {
+ int timestep = positionEvent->GetSender()->GetTimeStep();
+
+ mitk::Point3D click = positionEvent->GetWorldPosition();
+
+ mitk::Point3D first = this->m_Contour->GetVertexAt(0, timestep)->Coordinates;
+
+
+ if (first.EuclideanDistanceTo(click) < 1.5)
+ {
+ //finish
+ newStateEvent = new mitk::StateEvent(EIDYES, stateEvent->GetEvent());
+ }else
+ {
+ //stay active
+ newStateEvent = new mitk::StateEvent(EIDNO, stateEvent->GetEvent());
+ }
+ }
+
+ this->HandleEvent( newStateEvent );
+
+
+ return true;
+}
+
+
+
+
+bool mitk::LiveWireTool2D::OnFinish( Action* action, const StateEvent* stateEvent)
+{
+ // finish livewire tool interaction
+
+ /* check if event can be handled */
+ if ( Superclass::CanHandleEvent(stateEvent) < 1.0 ) return false;
+
+ const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
+ if (!positionEvent) return false;
+ /* END check if event can be handled */
+
+
+ //actual timestep
+ int timestep = positionEvent->GetSender()->GetTimeStep();
+
+ //remove last control point being added by double click
+ m_Contour->RemoveVertexAt(m_Contour->GetNumberOfVertices(timestep) - 1, timestep);
+
+ //save contour and corresponding plane geometry to list
+ std::pair contourPair(m_ContourModelNode.GetPointer(), dynamic_cast(positionEvent->GetSender()->GetCurrentWorldGeometry2D()->Clone().GetPointer()) );
+ m_Contours.push_back(contourPair);
+
+ m_LiveWireFilter->SetUseDynamicCostMap(false);
+ this->FinishTool();
+
+ return true;
+}
+
+
+
+void mitk::LiveWireTool2D::FinishTool()
+{
+
+ unsigned int numberOfTimesteps = m_Contour->GetTimeSlicedGeometry()->GetTimeSteps();
+
+ //close contour in each timestep
+ for( int i = 0; i <= numberOfTimesteps; i++)
+ {
+ m_Contour->Close(i);
+ }
+
+ //clear LiveWire node
+ m_ToolManager->GetDataStorage()->Remove( m_LiveWireContourNode );
+ m_LiveWireContourNode = NULL;
+ m_LiveWireContour = NULL;
+
+
+ //TODO visual feedback for completing livewire tool
+ m_ContourModelNode->AddProperty( "contour.color", ColorProperty::New(1.0, 1.0, 0.1), NULL, true );
+ m_ContourModelNode->SetProperty("name", StringProperty::New("contour node"));
+
+ //set the livewire interactor to edit control points
+ mitk::ContourModelLiveWireInteractor::Pointer interactor = mitk::ContourModelLiveWireInteractor::New(m_ContourModelNode);
+ interactor->SetWorkingImage(this->m_WorkingSlice);
+ m_ContourModelNode->SetInteractor(interactor);
+
+ //add interactor to globalInteraction instance
+ mitk::GlobalInteraction::GetInstance()->AddInteractor(interactor);
+ /* END complete livewire tool interaction */
+
+
+ /* reset contours and datanodes */
+ m_Contour = mitk::ContourModel::New();
+ m_ContourModelNode = mitk::DataNode::New();
+ m_ContourModelNode->SetData( m_Contour );
+ m_ContourModelNode->SetProperty("name", StringProperty::New("working contour node"));
+ m_ContourModelNode->SetProperty("visible", BoolProperty::New(true));
+ m_ContourModelNode->AddProperty( "contour.color", ColorProperty::New(0.9, 1.0, 0.1), NULL, true );
+ m_ContourModelNode->AddProperty( "points.color", ColorProperty::New(1.0, 0.0, 0.1), NULL, true );
+
+ m_LiveWireContour = mitk::ContourModel::New();
+ m_LiveWireContourNode = mitk::DataNode::New();
+ //m_LiveWireContourNode->SetData( m_LiveWireContour );
+ m_LiveWireContourNode->SetProperty("name", StringProperty::New("active livewire node"));
+ m_LiveWireContourNode->SetProperty("visible", BoolProperty::New(true));
+ m_LiveWireContourNode->AddProperty( "contour.color", ColorProperty::New(0.1, 1.0, 0.1), NULL, true );
+ m_LiveWireContourNode->AddProperty( "points.color", ColorProperty::New(0.5, 0.5, 0.1), NULL, true );
+ /* END reset contours and datanodes */
+}
+
+
+
+
+bool mitk::LiveWireTool2D::OnLastSegmentDelete( Action* action, const StateEvent* stateEvent)
+{
+ int timestep = stateEvent->GetEvent()->GetSender()->GetTimeStep();
+
+ //if last point of current contour will be removed go to start state and remove nodes
+ if( m_Contour->GetNumberOfVertices(timestep) <= 1 )
+ {
+ m_ToolManager->GetDataStorage()->Remove( m_LiveWireContourNode );
+ m_ToolManager->GetDataStorage()->Remove( m_ContourModelNode );
+ m_LiveWireContour = mitk::ContourModel::New();
+ m_Contour = mitk::ContourModel::New();
+ m_ContourModelNode->SetData( m_Contour );
+ m_LiveWireContourNode->SetData( m_LiveWireContour );
+ Superclass::Deactivated(); //go to start state
+ }
+ else //remove last segment from contour and reset livewire contour
+ {
+
+ m_LiveWireContour = mitk::ContourModel::New();
+
+ m_LiveWireContourNode->SetData(m_LiveWireContour);
+
+
+ mitk::ContourModel::Pointer newContour = mitk::ContourModel::New();
+ newContour->Expand(m_Contour->GetTimeSteps());
+
+ mitk::ContourModel::VertexIterator begin = m_Contour->IteratorBegin();
+
+ //iterate from last point to next active point
+ mitk::ContourModel::VertexIterator newLast = m_Contour->IteratorBegin() + (m_Contour->GetNumberOfVertices() - 1);
+
+ //go at least one down
+ if(newLast != begin)
+ {
+ newLast--;
+ }
+
+ //search next active control point
+ while(newLast != begin && !((*newLast)->IsControlPoint) )
+ {
+ newLast--;
+ }
+
+ //set position of start point for livewire filter to coordinates of the new last point
+ m_LiveWireFilter->SetStartPoint((*newLast)->Coordinates);
+
+ mitk::ContourModel::VertexIterator it = m_Contour->IteratorBegin();
+
+ //fill new Contour
+ while(it <= newLast)
+ {
+ newContour->AddVertex((*it)->Coordinates, (*it)->IsControlPoint, timestep);
+ it++;
+ }
+
+ newContour->SetIsClosed(m_Contour->IsClosed());
+
+ //set new contour visible
+ m_ContourModelNode->SetData(newContour);
+
+ m_Contour = newContour;
+
+ assert( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+ mitk::RenderingManager::GetInstance()->RequestUpdate( stateEvent->GetEvent()->GetSender()->GetRenderWindow() );
+ }
+
+ return true;
+}
+
+
+template
+void mitk::LiveWireTool2D::FindHighestGradientMagnitudeByITK(itk::Image* inputImage, itk::Index<3> &index, itk::Index<3> &returnIndex)
+{
+ typedef itk::Image InputImageType;
+ typedef typename InputImageType::IndexType IndexType;
+
+ unsigned long xMAX = inputImage->GetLargestPossibleRegion().GetSize()[0];
+ unsigned long yMAX = inputImage->GetLargestPossibleRegion().GetSize()[1];
+
+ returnIndex[0] = index[0];
+ returnIndex[1] = index[1];
+ returnIndex[2] = 0.0;
+
+
+ double gradientMagnitude = 0.0;
+ double maxGradientMagnitude = 0.0;
+
+ /*
+ the size and thus the region of 7x7 is only used to calculate the gradient magnitude in that region
+ not for searching the maximum value
+ */
+
+ //maximum value in each direction for size
+ typename InputImageType::SizeType size;
+ size[0] = 7;
+ size[1] = 7;
+
+ //minimum value in each direction for startRegion
+ IndexType startRegion;
+ startRegion[0] = index[0] - 3;
+ startRegion[1] = index[1] - 3;
+ if(startRegion[0] < 0) startRegion[0] = 0;
+ if(startRegion[1] < 0) startRegion[1] = 0;
+ if(xMAX - index[0] < 7) startRegion[0] = xMAX - 7;
+ if(yMAX - index[1] < 7) startRegion[1] = yMAX - 7;
+
+ index[0] = startRegion[0] + 3;
+ index[1] = startRegion[1] + 3;
+
+
+
+ typename InputImageType::RegionType region;
+ region.SetSize( size );
+ region.SetIndex( startRegion );
+
+ typedef typename itk::GradientMagnitudeImageFilter< InputImageType, InputImageType> GradientMagnitudeFilterType;
+ typename GradientMagnitudeFilterType::Pointer gradientFilter = GradientMagnitudeFilterType::New();
+ gradientFilter->SetInput(inputImage);
+ gradientFilter->GetOutput()->SetRequestedRegion(region);
+
+ gradientFilter->Update();
+ typename InputImageType::Pointer gradientMagnImage;
+ gradientMagnImage = gradientFilter->GetOutput();
+
+
+ IndexType currentIndex;
+ currentIndex[0] = 0;
+ currentIndex[1] = 0;
+
+ // search max (approximate) gradient magnitude
+ for( int x = -1; x <= 1; ++x)
+ {
+ currentIndex[0] = index[0] + x;
+
+ for( int y = -1; y <= 1; ++y)
+ {
+ currentIndex[1] = index[1] + y;
+
+ gradientMagnitude = gradientMagnImage->GetPixel(currentIndex);
+
+ //check for new max
+ if(maxGradientMagnitude < gradientMagnitude)
+ {
+ maxGradientMagnitude = gradientMagnitude;
+ returnIndex[0] = currentIndex[0];
+ returnIndex[1] = currentIndex[1];
+ returnIndex[2] = 0.0;
+ }//end if
+ }//end for y
+
+ currentIndex[1] = index[1];
+ }//end for x
+
+}
diff --git a/Modules/Segmentation/Interactions/mitkLiveWireTool2D.h b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.h
new file mode 100644
index 0000000000..99366a5b80
--- /dev/null
+++ b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.h
@@ -0,0 +1,121 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef mitkCorrectorTool2D_h_Included
+#define mitkCorrectorTool2D_h_Included
+
+#include "mitkCommon.h"
+#include "SegmentationExports.h"
+#include "mitkSegTool2D.h"
+#include
+#include
+
+#include
+
+namespace mitk
+{
+
+
+/**
+ \brief A 2D segmentation tool based on LiveWire approach.
+ The contour between the last user added point and the current mouse position
+ is computed by searching the shortest path according to specific features of
+ the image. The contour thus snappest to the boundary of objects.
+
+
+ \sa SegTool2D
+ \sa ImageLiveWireContourModelFilter
+
+ \ingroup Interaction
+ \ingroup ToolManagerEtAl
+
+ \warning Only to be instantiated by mitk::ToolManager.
+*/
+class Segmentation_EXPORT LiveWireTool2D : public SegTool2D
+{
+ public:
+
+ mitkClassMacro(LiveWireTool2D, SegTool2D);
+ itkNewMacro(LiveWireTool2D);
+
+ virtual const char** GetXPM() const;
+ virtual const char* GetName() const;
+
+ protected:
+
+ LiveWireTool2D();
+ virtual ~LiveWireTool2D();
+
+ /**
+ * \brief Calculates how good the data, this statemachine handles, is hit by the event.
+ *
+ */
+ virtual float CanHandleEvent( StateEvent const *stateEvent) const;
+
+ virtual void Activated();
+ virtual void Deactivated();
+
+ /// \brief Initialize tool
+ virtual bool OnInitLiveWire (Action*, const StateEvent*);
+
+ /// \brief Add a control point and finish current segment
+ virtual bool OnAddPoint (Action*, const StateEvent*);
+
+ /// \breif Actual LiveWire computation
+ virtual bool OnMouseMoved(Action*, const StateEvent*);
+
+ /// \brief Check double click on first control point to finish the LiveWire tool
+ virtual bool OnCheckPoint(Action*, const StateEvent*);
+
+ /// \brief Finish LiveWire tool
+ virtual bool OnFinish(Action*, const StateEvent*);
+
+ /// \brief Close the contour
+ virtual bool OnLastSegmentDelete(Action*, const StateEvent*);
+
+ /// \brief Don't use dynamic cost map for LiveWire calculation
+ virtual bool OnMouseMoveNoDynamicCosts(Action*, const StateEvent*);
+
+ /// \brief Finish contour interaction.
+ void FinishTool();
+
+
+ //the contour already set by the user
+ mitk::ContourModel::Pointer m_Contour;
+ //the corresponding datanode
+ mitk::DataNode::Pointer m_ContourModelNode;
+
+ //the current LiveWire computed contour
+ mitk::ContourModel::Pointer m_LiveWireContour;
+ //the corresponding datanode
+ mitk::DataNode::Pointer m_LiveWireContourNode;
+
+ //the current reference image
+ mitk::Image::Pointer m_WorkingSlice;
+
+ mitk::ImageLiveWireContourModelFilter::Pointer m_LiveWireFilter;
+
+ bool m_CreateAndUseDynamicCosts;
+
+ std::vector< std::pair > m_Contours;
+
+ template
+ void FindHighestGradientMagnitudeByITK(itk::Image* inputImage, itk::Index<3> &index, itk::Index<3> &returnIndex);
+};
+
+} // namespace
+
+#endif
diff --git a/Modules/Segmentation/Interactions/mitkLiveWireTool2D.xpm b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.xpm
new file mode 100644
index 0000000000..e5c4a22a1a
--- /dev/null
+++ b/Modules/Segmentation/Interactions/mitkLiveWireTool2D.xpm
@@ -0,0 +1,99 @@
+/* XPM */
+static const char * mitkLiveWireTool2D_xpm[] = {
+"24 24 72 1",
+" c None",
+". c #179FBF",
+"+ c #01AAD2",
+"@ c #00AAD3",
+"# c #0BA4CA",
+"$ c #14A0C1",
+"% c #11A2C4",
+"& c #06A7CE",
+"* c #0AA5CB",
+"= c #2895B1",
+"- c #229AB5",
+"; c #02A9D1",
+"> c #0DA3C8",
+", c #1F99B7",
+"' c #2598B5",
+") c #0CA4C9",
+"! c #13A1C2",
+"~ c #03A8D0",
+"{ c #2795B1",
+"] c #179EBF",
+"^ c #04A9D0",
+"/ c #1C9BBA",
+"( c #2497B3",
+"_ c #1E9AB9",
+": c #209AB9",
+"< c #13A0C3",
+"[ c #03A9D1",
+"} c #189DBE",
+"| c #07A7CD",
+"1 c #1A9CBC",
+"2 c #1E9BBA",
+"3 c #2699B4",
+"4 c #2498B4",
+"5 c #229AB6",
+"6 c #0EA3C8",
+"7 c #07A7CE",
+"8 c #189EBE",
+"9 c #1C9CBA",
+"0 c #04A7CF",
+"a c #11A1C5",
+"b c #07A6CD",
+"c c #09A6CC",
+"d c #2A94AD",
+"e c #1C9CBB",
+"f c #1E9BB8",
+"g c #259AB4",
+"h c #05A7CE",
+"i c #2E92A9",
+"j c #0EA2C6",
+"k c #3A8FA2",
+"l c #1B9CBC",
+"m c #1D9BB9",
+"n c #169FC1",
+"o c #00A9D2",
+"p c #1D9BBA",
+"q c #2B94AC",
+"r c #16A0C1",
+"s c #10A2C6",
+"t c #209BB9",
+"u c #1D9CBB",
+"v c #179FC0",
+"w c #000000",
+"x c #219AB4",
+"y c #FFFFFF",
+"z c #F3F6F1",
+"A c #189EC0",
+"B c #00AAD2",
+"C c #16A7BC",
+"D c #17A0BF",
+"E c #2299B7",
+"F c #1D9BB8",
+"G c #1A9CBD",
+" .+@#$%&@*= ",
+" -;>, '); ",
+" !~{ &) ",
+" ]^ /@ ",
+" (@_ :@ ",
+" <[ }| ",
+" @1 +2 ",
+" @345 67 ",
+" +@890 a^ ",
+" }bcdef gh|i ",
+" jk&@l(mn^op ",
+" qr+@@@@st ",
+" u& ",
+" 0vwwwwwwwww ",
+" @xwwwyyzyzww ",
+" ABwwwwzzzzzzww ",
+" oCwwwwwwwwzzzw ",
+" D@ wwzzw ",
+" c@@@ wzzw ",
+" EFG wwzzw ",
+" wwwwwwwwzzzw ",
+" wwwwzzzzzyww ",
+" wwwwyyyyyww ",
+" wwwwwwwwww "};
diff --git a/Modules/Segmentation/Interactions/mitkSegTool2D.cpp b/Modules/Segmentation/Interactions/mitkSegTool2D.cpp
index 739ed49b80..e4abb59b68 100644
--- a/Modules/Segmentation/Interactions/mitkSegTool2D.cpp
+++ b/Modules/Segmentation/Interactions/mitkSegTool2D.cpp
@@ -1,373 +1,414 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "mitkSegTool2D.h"
#include "mitkToolManager.h"
#include "mitkDataStorage.h"
#include "mitkBaseRenderer.h"
#include "mitkPlaneGeometry.h"
#include "mitkExtractImageFilter.h"
#include "mitkExtractDirectedPlaneImageFilter.h"
//Include of the new ImageExtractor
#include "mitkExtractDirectedPlaneImageFilterNew.h"
#include "mitkPlanarCircle.h"
#include "mitkOverwriteSliceImageFilter.h"
#include "mitkOverwriteDirectedPlaneImageFilter.h"
#include "mitkGetModuleContext.h"
//Includes for 3DSurfaceInterpolation
#include "mitkImageToContourFilter.h"
#include "mitkSurfaceInterpolationController.h"
//includes for resling and overwriting
#include
#include
#include
#include
#include
#include "mitkOperationEvent.h"
#include "mitkUndoController.h"
#define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a)))
mitk::SegTool2D::SegTool2D(const char* type)
:Tool(type),
- m_LastEventSender(NULL),
- m_LastEventSlice(0),
- m_Contourmarkername ("Position"),
- m_ShowMarkerNodes (false),
- m_3DInterpolationEnabled(true)
+m_LastEventSender(NULL),
+m_LastEventSlice(0),
+m_Contourmarkername ("Position"),
+m_ShowMarkerNodes (false),
+m_3DInterpolationEnabled(true)
{
}
mitk::SegTool2D::~SegTool2D()
{
}
float mitk::SegTool2D::CanHandleEvent( StateEvent const *stateEvent) const
{
const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent());
if (!positionEvent) return 0.0;
if ( positionEvent->GetSender()->GetMapperID() != BaseRenderer::Standard2D ) return 0.0; // we don't want anything but 2D
//This are the mouse event that are used by the statemachine patterns for zooming and panning. This must be possible although a tool is activ
if (stateEvent->GetId() == EIDRIGHTMOUSEBTN || stateEvent->GetId() == EIDMIDDLEMOUSEBTN || stateEvent->GetId() == EIDRIGHTMOUSEBTNANDCTRL ||
- stateEvent->GetId() == EIDMIDDLEMOUSERELEASE || stateEvent->GetId() == EIDRIGHTMOUSERELEASE || stateEvent->GetId() == EIDRIGHTMOUSEBTNANDMOUSEMOVE ||
- stateEvent->GetId() == EIDMIDDLEMOUSEBTNANDMOUSEMOVE || stateEvent->GetId() == EIDCTRLANDRIGHTMOUSEBTNANDMOUSEMOVE || stateEvent->GetId() == EIDCTRLANDRIGHTMOUSEBTNRELEASE )
+ stateEvent->GetId() == EIDMIDDLEMOUSERELEASE || stateEvent->GetId() == EIDRIGHTMOUSERELEASE || stateEvent->GetId() == EIDRIGHTMOUSEBTNANDMOUSEMOVE ||
+ stateEvent->GetId() == EIDMIDDLEMOUSEBTNANDMOUSEMOVE || stateEvent->GetId() == EIDCTRLANDRIGHTMOUSEBTNANDMOUSEMOVE || stateEvent->GetId() == EIDCTRLANDRIGHTMOUSEBTNRELEASE )
{
//Since the usual segmentation tools currently do not need right click interaction but the mitkDisplayVectorInteractor
return 0.0;
}
else
{
return 1.0;
}
}
bool mitk::SegTool2D::DetermineAffectedImageSlice( const Image* image, const PlaneGeometry* plane, int& affectedDimension, int& affectedSlice )
{
assert(image);
assert(plane);
// compare normal of plane to the three axis vectors of the image
Vector3D normal = plane->GetNormal();
Vector3D imageNormal0 = image->GetSlicedGeometry()->GetAxisVector(0);
Vector3D imageNormal1 = image->GetSlicedGeometry()->GetAxisVector(1);
Vector3D imageNormal2 = image->GetSlicedGeometry()->GetAxisVector(2);
normal.Normalize();
imageNormal0.Normalize();
imageNormal1.Normalize();
imageNormal2.Normalize();
imageNormal0.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal0.Get_vnl_vector()) );
imageNormal1.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal1.Get_vnl_vector()) );
imageNormal2.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal2.Get_vnl_vector()) );
double eps( 0.00001 );
// axial
if ( imageNormal2.GetNorm() <= eps )
{
affectedDimension = 2;
}
// sagittal
else if ( imageNormal1.GetNorm() <= eps )
{
affectedDimension = 1;
}
// frontal
else if ( imageNormal0.GetNorm() <= eps )
{
affectedDimension = 0;
}
else
{
affectedDimension = -1; // no idea
return false;
}
// determine slice number in image
Geometry3D* imageGeometry = image->GetGeometry(0);
Point3D testPoint = imageGeometry->GetCenter();
Point3D projectedPoint;
plane->Project( testPoint, projectedPoint );
Point3D indexPoint;
imageGeometry->WorldToIndex( projectedPoint, indexPoint );
affectedSlice = ROUND( indexPoint[affectedDimension] );
MITK_DEBUG << "indexPoint " << indexPoint << " affectedDimension " << affectedDimension << " affectedSlice " << affectedSlice;
// check if this index is still within the image
if ( affectedSlice < 0 || affectedSlice >= static_cast(image->GetDimension(affectedDimension)) ) return false;
return true;
}
+
mitk::Image::Pointer mitk::SegTool2D::GetAffectedImageSliceAs2DImage(const PositionEvent* positionEvent, const Image* image)
{
if (!positionEvent) return NULL;
assert( positionEvent->GetSender() ); // sure, right?
unsigned int timeStep = positionEvent->GetSender()->GetTimeStep( image ); // get the timestep of the visible part (time-wise) of the image
// first, we determine, which slice is affected
const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) );
+ return this->GetAffectedImageSliceAs2DImage(planeGeometry, image, timeStep);
+}
+
+
+mitk::Image::Pointer mitk::SegTool2D::GetAffectedImageSliceAs2DImage(const PlaneGeometry* planeGeometry, const Image* image, unsigned int timeStep)
+{
if ( !image || !planeGeometry ) return NULL;
//Make sure that for reslicing and overwriting the same alogrithm is used. We can specify the mode of the vtk reslicer
vtkSmartPointer reslice = vtkSmartPointer