From 81e5245951244bd6df71c1069f24516aca486f61 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Tue, 18 Oct 2016 18:01:11 -0500 Subject: [PATCH] ... --- lib/fpa/Image/SkeletonFilter.h | 10 ++- lib/fpa/Image/SkeletonFilter.hxx | 113 ++++++++++++++------------- plugins/Plugins/SkeletonFilter.cxx | 118 +++++++++++++++++++++++++++++ plugins/Plugins/SkeletonFilter.h | 32 ++++++++ 4 files changed, 219 insertions(+), 54 deletions(-) create mode 100644 plugins/Plugins/SkeletonFilter.cxx create mode 100644 plugins/Plugins/SkeletonFilter.h diff --git a/lib/fpa/Image/SkeletonFilter.h b/lib/fpa/Image/SkeletonFilter.h index 4d36fd1..5951709 100644 --- a/lib/fpa/Image/SkeletonFilter.h +++ b/lib/fpa/Image/SkeletonFilter.h @@ -35,6 +35,7 @@ namespace fpa typedef cpExtensions::DataStructures::Skeleton< TCostMap::ImageDimension > TSkeleton; + typedef itk::Image< unsigned char, _TCostMap::ImageDimension > TMarks; typedef itk::Functor::IndexLexicographicCompare< _TCostMap::ImageDimension > @@ -68,6 +69,9 @@ namespace fpa TSkeleton* GetSkeleton( ); const TSkeleton* GetSkeleton( ) const; + TMarks* GetMarks( ); + const TMarks* GetMarks( ) const; + protected: SkeletonFilter( ); virtual ~SkeletonFilter( ); @@ -78,15 +82,15 @@ namespace fpa const TDistanceMap* dmap, const TCostMap* cmap, const TMST* mst, - TIndices* end_points + TIndicesData& end_points ); void _Skeleton( const TDistanceMap* dmap, const TCostMap* cmap, const TMST* mst, - TIndices* end_points - TIndices* bifurcations, + const TIndicesData& end_points, + TIndicesData& bifurcations, TSkeleton* skeleton ); diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index 65a3fe2..9488ed8 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -1,72 +1,67 @@ #ifndef __fpa__Image__SkeletonFilter__hxx__ #define __fpa__Image__SkeletonFilter__hxx__ +#include +#include + // ------------------------------------------------------------------------- -#define __fpa__Image__SkeletonFilter__InputMacro( input_name, input_type, input_id ) \ +#define fpa_Image_SkeletonFilter_InputMacro( i_n, i_t, i_i ) \ template< class _TDistanceMap, class _TCostMap > \ typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - input_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##input_name( ) \ + i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##i_n( ) \ { \ return( \ - dynamic_cast< input_type* >( \ - this->Superclass::GetInput( input_id ) \ - ) \ + dynamic_cast< i_t* >( this->Superclass::GetInput( i_i ) ) \ ); \ } \ template< class _TDistanceMap, class _TCostMap > \ const \ typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - input_type* \ - fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##input_name( ) const \ + i_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##i_n( ) const \ { \ return( \ - dynamic_cast< const input_type* >( \ - this->Superclass::GetInput( input_id ) \ - ) \ + dynamic_cast< const i_t* >( this->Superclass::GetInput( i_i ) ) \ ); \ } \ template< class _TDistanceMap, class _TCostMap > \ void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Set##input_name( input_type* input ) \ + Set##i_n( i_t* input ) \ { \ - this->Superclass::SetInput( input_id, input ); \ + this->Superclass::SetNthInput( i_i, input ); \ } -__fpa__Image__SkeletonFilter__InputMacro( DistanceMap, TDistanceMap, 0 ); -__fpa__Image__SkeletonFilter__InputMacro( CostMap, TCostMap, 1 ); -__fpa__Image__SkeletonFilter__InputMacro( MinimumSpanningTree, TMST, 2 ); +fpa_Image_SkeletonFilter_InputMacro( DistanceMap, TDistanceMap, 0 ); +fpa_Image_SkeletonFilter_InputMacro( CostMap, TCostMap, 1 ); +fpa_Image_SkeletonFilter_InputMacro( MinimumSpanningTree, TMST, 2 ); // ------------------------------------------------------------------------- -#define __fpa__Image__SkeletonFilter__OutputMacro( output_name, output_type, output_id ) \ +#define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t, o_i ) \ template< class _TDistanceMap, class _TCostMap > \ typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##output_name( ) \ + o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##o_n( ) \ { \ return( \ - dynamic_cast< output_type* >( \ - this->Superclass::GetOutput( output_id ) \ - ) \ + dynamic_cast< o_t* >( this->Superclass::GetOutput( o_i ) ) \ ); \ } \ template< class _TDistanceMap, class _TCostMap > \ const typename \ fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ - Get##output_name( ) const \ + o_t* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##o_n( ) const \ { \ return( \ - dynamic_cast< const output_type* >( \ - this->Superclass::GetOutput( output_id ) \ - ) \ + dynamic_cast< const o_t* >( this->Superclass::GetOutput( o_i ) ) \ ); \ } -__fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 ); -__fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 ); -__fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 ); +fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton, 0 ); +fpa_Image_SkeletonFilter_OutputMacro( EndPoints, TIndices, 1 ); +fpa_Image_SkeletonFilter_OutputMacro( Bifurcations, TIndices, 2 ); +fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks, 3 ); // ------------------------------------------------------------------------- template< class _TDistanceMap, class _TCostMap > @@ -74,15 +69,17 @@ fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: SkeletonFilter( ) : Superclass( ) { - this->SetNumbeOfRequiredInputs( 3 ); - this->SetNumbeOfRequiredOutputs( 3 ); + this->SetNumberOfRequiredInputs( 3 ); + this->SetNumberOfRequiredOutputs( 3 ); typename TIndices::Pointer end_points = TIndices::New( ); typename TIndices::Pointer bifurcations = TIndices::New( ); typename TSkeleton::Pointer skeleton = TSkeleton::New( ); - this->SetNthOutput( 0, skeleton ); - this->SetNthOutput( 1, end_points ); - this->SetNthOutput( 2, bifurcations ); + typename TMarks::Pointer marks = TMarks::New( ); + this->SetNthOutput( 0, skeleton.GetPointer( ) ); + this->SetNthOutput( 1, end_points.GetPointer( ) ); + this->SetNthOutput( 2, bifurcations.GetPointer( ) ); + this->SetNthOutput( 3, marks.GetPointer( ) ); } // ------------------------------------------------------------------------- @@ -101,18 +98,18 @@ GenerateData( ) const TDistanceMap* dmap = this->GetDistanceMap( ); const TCostMap* cmap = this->GetCostMap( ); const TMST* mst = this->GetMinimumSpanningTree( ); - TIndices* end_points = this->GetEndPoints( ); - TIndices* bifurcations = this->GetBifurcations( ); - TSkeleton* skeleton = this->GetSkeleton( ); + TIndices* ep = this->GetEndPoints( ); + TIndices* bi = this->GetBifurcations( ); + TSkeleton* sk = this->GetSkeleton( ); // 1. Check input correspondance // TODO // 2. Detect end-points - this->_EndPoints( dmap, cmap, mst, end_points ); + this->_EndPoints( dmap, cmap, mst, ep->Get( ) ); // 3. Build skeleton and keep track of bifurcations - this->_Skeleton( dmap, cmap, mst, end_points, bifurcations, skeleton ); + this->_Skeleton( dmap, cmap, mst, ep->Get( ), bi->Get( ), sk ); } // ------------------------------------------------------------------------- @@ -122,7 +119,7 @@ _EndPoints( const TDistanceMap* dmap, const TCostMap* cmap, const TMST* mst, - TIndices* end_points + TIndicesData& end_points ) { typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt; @@ -132,7 +129,8 @@ _EndPoints( typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks; // Some values - typename _TMarks::Pointer marks = _TMarks::New( ); + // typename _TMarks::Pointer marks = _TMarks::New( ); + auto marks = this->GetMarks( ); marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) ); marks->SetRequestedRegion( dmap->GetRequestedRegion( ) ); marks->SetBufferedRegion( dmap->GetBufferedRegion( ) ); @@ -161,6 +159,8 @@ _EndPoints( } // rof // BFS from maximum queue + auto region = dmap->GetRequestedRegion( ); + double init_v = queue.begin( )->first; while( queue.size( ) > 0 ) { // Get node @@ -169,24 +169,26 @@ _EndPoints( queue.erase( nIt ); unsigned char m = marks->GetPixel( n.second ); - if( ( m & 0x01 ) == 0x01 ) + // if( ( m & 0x01 ) == 0x01 ) + if( m != 0 || ( n.first / init_v ) < double( 1e-1 ) ) continue; + std::cout << n.second << " " << n.first << " " << dmap->GetPixel( n.second ) << std::endl; + // Mark it and update end-points m |= 0x01; marks->SetPixel( n.second, m ); - end_points->Get( ).insert( n.second ); + end_points.insert( n.second ); // Get path typename TMST::TPath::Pointer path; mst->GetPath( path, n.second ); for( unsigned long i = 0; i < path->GetSize( ); ++i ) { - TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i ); + TIndex idx = path->GetVertex( i ); typename _TCostMap::PointType cnt; cmap->TransformIndexToPhysicalPoint( idx, cnt ); - double d = double( dmap->GetPixel( idx ) ); - d *= double( 2 ); + double d = double( 2 ) * double( dmap->GetPixel( idx ) ); // Mark sphere std::queue< TIndex > q; @@ -198,7 +200,7 @@ _EndPoints( unsigned char m = marks->GetPixel( v ); if( ( m & 0x02 ) == 0x02 ) continue; - m |= 0x03; + m |= 0x02; marks->SetPixel( v, m ); for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x ) @@ -234,11 +236,20 @@ _Skeleton( const TDistanceMap* dmap, const TCostMap* cmap, const TMST* mst, - TIndices* end_points - TIndices* bifurcations, + const TIndicesData& end_points, + TIndicesData& bifurcations, TSkeleton* skeleton ) { + typedef typename TMST::TPath _TPath; + + for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) + { + typename _TPath::Pointer path; + mst->GetPath( path, *eIt ); + skeleton->AddBranch( path ); + + } // rof } #endif // __fpa__Image__SkeletonFilter__hxx__ diff --git a/plugins/Plugins/SkeletonFilter.cxx b/plugins/Plugins/SkeletonFilter.cxx new file mode 100644 index 0000000..f00a0bf --- /dev/null +++ b/plugins/Plugins/SkeletonFilter.cxx @@ -0,0 +1,118 @@ +#include +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +fpaPlugins::SkeletonFilter:: +SkeletonFilter( ) + : Superclass( ) +{ + typedef cpPlugins::DataObjects::Image _TImage; + typedef cpPlugins::DataObjects::Skeleton _TSkeleton; + + this->_ConfigureInput< _TImage >( "DistanceMap", true, false ); + this->_ConfigureInput< _TImage >( "CostMap", true, false ); + this->_ConfigureInput< _TImage >( "MST", true, false ); + this->_ConfigureOutput< _TSkeleton >( "Skeleton" ); + this->_ConfigureOutput< _TImage >( "Marks" ); + /* TODO + this->_ConfigureOutput< _TMesh >( "EndPoints" ); + this->_ConfigureOutput< _TMesh >( "Bifurcations" ); + */ +} + +// ------------------------------------------------------------------------- +fpaPlugins::SkeletonFilter:: +~SkeletonFilter( ) +{ +} + +// ------------------------------------------------------------------------- +void fpaPlugins::SkeletonFilter:: +_GenerateData( ) +{ + auto o = this->GetInputData( "DistanceMap" ); + cpPlugins_Demangle_ImageScalars_Dims( o, _GD0 ); + else this->_Error( "Invalid input image." ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap > +void fpaPlugins::SkeletonFilter:: +_GD0( _TDistanceMap* dmap ) +{ + auto cmap = this->GetInputData< _TDistanceMap >( "CostMap" ); + if( cmap != NULL ) + this->_GD1( dmap, cmap ); + else + this->_Error( "Temporary error: invalid cost map." ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +void fpaPlugins::SkeletonFilter:: +_GD1( _TDistanceMap* dmap, _TCostMap* cmap ) +{ + typedef fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap > _TFilter; + typedef typename _TFilter::TMST _TMST; + + auto mst = this->GetInputData< _TMST >( "MST" ); + if( mst == NULL ) + this->_Error( "Invalid MST." ); + + auto filter = this->_CreateITK< _TFilter >( ); + filter->SetDistanceMap( dmap ); + filter->SetCostMap( cmap ); + filter->SetMinimumSpanningTree( mst ); + filter->Update( ); + this->GetOutput( "Skeleton" )->SetITK( filter->GetSkeleton( ) ); + this->GetOutput( "Marks" )->SetITK( filter->GetMarks( ) ); + + + /* TODO + auto ep = filter->GetEndPoints( ); + auto bi = filter->GetBifurcations( ); + + auto ep_pd = this->GetOutputData< vtkPolyData >( "EndPoints" ); + if( ep_pd == NULL ) + { + auto points = vtkSmartPointer< vtkPoints >::New( ); + auto verts = vtkSmartPointer< vtkCellArray >::New( ); + auto lines = vtkSmartPointer< vtkCellArray >::New( ); + auto polys = vtkSmartPointer< vtkCellArray >::New( ); + auto strips = vtkSmartPointer< vtkCellArray >::New( ); + auto pd = vtkSmartPointer< vtkPolyData >::New( ); + pd->SetPoints( points ); + pd->SetVerts( verts ); + pd->SetLines( lines ); + pd->SetPolys( polys ); + pd->SetStrips( strips ); + + this->GetOutput( "EndPoints" )->SetVTK( pd ); + ep_pd = this->GetOutputData< vtkPolyData >( "EndPoints" ); + + } // fi + + for( auto iIt = ep.begin( ); iIt != ep.end( ); ++iIt ) + { + typename _TCostMap::PointType p; + cmap->TransformIndexToPhysicalPoint( *iIt, p ); + + if( _TCostMap::ImageDimension == 1 ) + ep_pd->GetPoints( )->InsertNextPoint( p[ 0 ], 0, 0 ); + else if( _TCostMap::ImageDimension == 2 ) + ep_pd->GetPoints( )->InsertNextPoint( p[ 0 ], p[ 1 ], 0 ); + else if( _TCostMap::ImageDimension > 2 ) + ep_pd->GetPoints( )->InsertNextPoint( p[ 0 ], p[ 1 ], p[ 2 ] ); + + ep_pd->GetVerts( )->InsertNextCell( 1 ); + ep_pd->GetVerts( )->InsertCellPoint( ep_pd->GetNumberOfPoints( ) - 1 ); + + } // rof + */ +} + +// eof - $RCSfile$ diff --git a/plugins/Plugins/SkeletonFilter.h b/plugins/Plugins/SkeletonFilter.h new file mode 100644 index 0000000..a24de80 --- /dev/null +++ b/plugins/Plugins/SkeletonFilter.h @@ -0,0 +1,32 @@ +#ifndef __fpa__Plugins__SkeletonFilter__h__ +#define __fpa__Plugins__SkeletonFilter__h__ + +#include +#include + +namespace fpaPlugins +{ + /** + */ + class fpaPlugins_EXPORT SkeletonFilter + : public cpPlugins::BaseObjects::ProcessObject + { + cpPluginsObject( + SkeletonFilter, + cpPlugins::BaseObjects::ProcessObject, + fpa + ); + + protected: + template< class _TDistanceMap > + inline void _GD0( _TDistanceMap* dmap ); + + template< class _TDistanceMap, class _TCostMap > + inline void _GD1( _TDistanceMap* dmap, _TCostMap* cmap ); + }; + +} // ecapseman + +#endif // __fpa__Plugins__SkeletonFilter__h__ + +// eof - $RCSfile$ -- 2.45.1