From 9622bd5b833a8845881003228207e0caca59b081 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Tue, 9 Dec 2014 09:58:31 +0100 Subject: [PATCH] First commit --- CMakeLists.txt | 96 ++++++ COMPILATION | 27 ++ README | 21 ++ appli/CMakeLists.txt | 6 + appli/examples/CMakeLists.txt | 31 ++ .../example_ImageAlgorithmDijkstra_00.cxx | 44 +++ .../example_ImageAlgorithmDijkstra_01.cxx | 152 +++++++++ .../example_ImageAlgorithmDijkstra_02.cxx | 205 ++++++++++++ .../example_ImageAlgorithmFastMarching_00.cxx | 44 +++ .../example_ImageAlgorithmFastMarching_01.cxx | 152 +++++++++ .../example_ImageAlgorithmRegionGrow_00.cxx | 51 +++ .../example_ImageAlgorithmRegionGrow_01.cxx | 157 +++++++++ ...AlgorithmRegionGrow_MultipleThresholds.cxx | 118 +++++++ cmake/CMakeLists.txt | 7 + cmake/FrontAlgorithmsConfig.cmake.in | 43 +++ data/binary_test_2D_00.png | Bin 0 -> 5192 bytes data/ones_image.png | Bin 0 -> 1697 bytes lib/CMakeLists.txt | 67 ++++ lib/fpa/Base/Algorithm.h | 153 +++++++++ lib/fpa/Base/Algorithm.hxx | 300 ++++++++++++++++++ lib/fpa/Base/Dijkstra.h | 120 +++++++ lib/fpa/Base/Dijkstra.hxx | 95 ++++++ lib/fpa/Base/Events.h | 136 ++++++++ lib/fpa/Base/FastMarching.h | 76 +++++ lib/fpa/Base/FastMarching.hxx | 121 +++++++ lib/fpa/Base/Functors/InvertCostFunction.h | 56 ++++ lib/fpa/Base/RegionGrow.h | 117 +++++++ lib/fpa/Base/RegionGrow.hxx | 87 +++++ lib/fpa/Base/TreeExtractor.h | 84 +++++ lib/fpa/Base/TreeExtractor.hxx | 250 +++++++++++++++ lib/fpa/Common.cxx.in | 9 + lib/fpa/Common.h | 15 + lib/fpa/Image/Algorithm.h | 105 ++++++ lib/fpa/Image/Algorithm.hxx | 180 +++++++++++ lib/fpa/Image/Dijkstra.h | 58 ++++ lib/fpa/Image/FastMarching.h | 58 ++++ .../Functors/RegionGrowAllBelongsFunction.h | 78 +++++ .../Functors/RegionGrowThresholdFunction.h | 96 ++++++ lib/fpa/Image/RegionGrow.h | 83 +++++ .../Image/RegionGrowWithMultipleThresholds.h | 70 ++++ .../RegionGrowWithMultipleThresholds.hxx | 137 ++++++++ lib/fpa/VTK/Image2DObserver.h | 91 ++++++ lib/fpa/VTK/Image2DObserver.hxx | 188 +++++++++++ lib/fpa/VTK/Image3DObserver.h | 93 ++++++ lib/fpa/VTK/Image3DObserver.hxx | 203 ++++++++++++ lib/fpa/VTK/ImageMPR.cxx | 159 ++++++++++ lib/fpa/VTK/ImageMPR.h | 67 ++++ 47 files changed, 4506 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 COMPILATION create mode 100644 README create mode 100644 appli/CMakeLists.txt create mode 100644 appli/examples/CMakeLists.txt create mode 100644 appli/examples/example_ImageAlgorithmDijkstra_00.cxx create mode 100644 appli/examples/example_ImageAlgorithmDijkstra_01.cxx create mode 100644 appli/examples/example_ImageAlgorithmDijkstra_02.cxx create mode 100644 appli/examples/example_ImageAlgorithmFastMarching_00.cxx create mode 100644 appli/examples/example_ImageAlgorithmFastMarching_01.cxx create mode 100644 appli/examples/example_ImageAlgorithmRegionGrow_00.cxx create mode 100644 appli/examples/example_ImageAlgorithmRegionGrow_01.cxx create mode 100644 appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx create mode 100644 cmake/CMakeLists.txt create mode 100644 cmake/FrontAlgorithmsConfig.cmake.in create mode 100644 data/binary_test_2D_00.png create mode 100644 data/ones_image.png create mode 100644 lib/CMakeLists.txt create mode 100644 lib/fpa/Base/Algorithm.h create mode 100644 lib/fpa/Base/Algorithm.hxx create mode 100644 lib/fpa/Base/Dijkstra.h create mode 100644 lib/fpa/Base/Dijkstra.hxx create mode 100644 lib/fpa/Base/Events.h create mode 100644 lib/fpa/Base/FastMarching.h create mode 100644 lib/fpa/Base/FastMarching.hxx create mode 100644 lib/fpa/Base/Functors/InvertCostFunction.h create mode 100644 lib/fpa/Base/RegionGrow.h create mode 100644 lib/fpa/Base/RegionGrow.hxx create mode 100644 lib/fpa/Base/TreeExtractor.h create mode 100644 lib/fpa/Base/TreeExtractor.hxx create mode 100644 lib/fpa/Common.cxx.in create mode 100644 lib/fpa/Common.h create mode 100644 lib/fpa/Image/Algorithm.h create mode 100644 lib/fpa/Image/Algorithm.hxx create mode 100644 lib/fpa/Image/Dijkstra.h create mode 100644 lib/fpa/Image/FastMarching.h create mode 100644 lib/fpa/Image/Functors/RegionGrowAllBelongsFunction.h create mode 100644 lib/fpa/Image/Functors/RegionGrowThresholdFunction.h create mode 100644 lib/fpa/Image/RegionGrow.h create mode 100644 lib/fpa/Image/RegionGrowWithMultipleThresholds.h create mode 100644 lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx create mode 100644 lib/fpa/VTK/Image2DObserver.h create mode 100644 lib/fpa/VTK/Image2DObserver.hxx create mode 100644 lib/fpa/VTK/Image3DObserver.h create mode 100644 lib/fpa/VTK/Image3DObserver.hxx create mode 100644 lib/fpa/VTK/ImageMPR.cxx create mode 100644 lib/fpa/VTK/ImageMPR.h diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..cd777b0 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,96 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.6) + +# for CMake 2.6 corrected behaviour (see "cmake --help-policy CMP0003") +IF( + COMMAND cmake_policy AND + ${CMAKE_MAJOR_VERSION} EQUAL 2 AND + ${CMAKE_MINOR_VERSION} GREATER 4 + ) + CMAKE_POLICY(SET CMP0003 NEW) + CMAKE_POLICY(SET CMP0005 NEW) + CMAKE_POLICY(SET CMP0011 NEW) + CMAKE_POLICY(SET CMP0012 NEW) +ENDIF( + COMMAND cmake_policy AND + ${CMAKE_MAJOR_VERSION} EQUAL 2 AND + ${CMAKE_MINOR_VERSION} GREATER 4 + ) + +## ================ +## = Project name = +## ================ + +PROJECT(FrontAlgorithms) +SET(FrontAlgorithms_MAJOR_VERSION "0") +SET(FrontAlgorithms_MINOR_VERSION "0") +SET(FrontAlgorithms_RELEASE_VERSION "1") +SET(FrontAlgorithms_VERSION "${FrontAlgorithms_MAJOR_VERSION}.${FrontAlgorithms_MINOR_VERSION}.${FrontAlgorithms_RELEASE_VERSION}") + +## =========== +## = Options = +## =========== + +OPTION(BUILD_EXAMPLES "Build examples" OFF) +OPTION(BUILD_SHARED_LIBS "Build shared libs" OFF) +OPTION(USE_VTK "Build using VTK" OFF) + +IF(BUILD_SHARED_LIBS) + SET(LIB_TYPE SHARED) +ELSE(BUILD_SHARED_LIBS) + SET(LIB_TYPE STATIC) +ENDIF(BUILD_SHARED_LIBS) + +## ============ +## = Packages = +## ============ + +INCLUDE(GenerateExportHeader) + +FIND_PACKAGE(ITK REQUIRED) +INCLUDE(${ITK_USE_FILE}) + +IF(USE_VTK) + FIND_PACKAGE(VTK REQUIRED) + INCLUDE(${VTK_USE_FILE}) +ENDIF(USE_VTK) + +## ================================================ +## = Do not allow to build inside the source tree = +## ================================================ + +IF(PROJECT_BINARY_DIR STREQUAL ${PROJECT_SOURCE_DIR}) + MESSAGE(FATAL_ERROR "Building in the source tree is not allowed") +ENDIF(PROJECT_BINARY_DIR STREQUAL ${PROJECT_SOURCE_DIR}) + +## ===================================== +## = Where to put executables and libs = +## ===================================== + +SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}) +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}) +MARK_AS_ADVANCED( + CMAKE_BACKWARDS_COMPATIBILITY + EXECUTABLE_OUTPUT_PATH + LIBRARY_OUTPUT_PATH + ) + +## ============================== +## = Subdirs containing headers = +## ============================== + +INCLUDE_DIRECTORIES( + ${PROJECT_SOURCE_DIR}/lib + ${PROJECT_BINARY_DIR}/lib + ) + +## =========================== +## = Subdirs containing code = +## =========================== + +SUBDIRS( + cmake + lib + appli + ) + +## eof - $RCSfile$ diff --git a/COMPILATION b/COMPILATION new file mode 100644 index 0000000..1886dd8 --- /dev/null +++ b/COMPILATION @@ -0,0 +1,27 @@ +@description + The project uses CMake as compilation helper. It should compile on the three + major box flavors (linux, windows, mac). However, as of december 2014, it + has only been tested on linux Fedora 20 and Ubuntu 14.04. + +@prerequisites + 1. CMake (>=2.8.12.2) + + 3. Insight Toolkit -ITK- (>=4.6.0) + 3.1 Required cmake flags: + BUILD_SHARED_LIBS:BOOL=ON + 3.2 [OPTIONAL] If USE_VTK=ON + Module_ITKVtkGlue:BOOL=ON + + 4. [OPTIONAL] Visualization Toolkit -VTK- (>=6.1.0) + 4.1 Required cmake flags: + BUILD_SHARED_LIBS:BOOL=ON + +@cmake_flags + BUILD_EXAMPLES:BOOL + Build example applications? (most of them are command line) + CMAKE_BUILD_TYPE:STRING + Debug/Release? + CMAKE_INSTALL_PREFIX:STRING + Where to put installation products? (in windows this option has no use) + +## eof - $RCSfile$ diff --git a/README b/README new file mode 100644 index 0000000..4ebdcfb --- /dev/null +++ b/README @@ -0,0 +1,21 @@ + +@project + FrontAlgorithms: Generic implementation of front propagation algorithms + with some extra features + +@version + 0.0.1 (2014-12-31) + +@authors + Maciej ORKISZ (maciej.orkisz@creatis.insa-lyon.fr) + Leonardo FLOREZ-VALENCIA (florez-l@javeriana.edu.co) + + +@description + + +@license + + + +## eof - $RCSfile$ diff --git a/appli/CMakeLists.txt b/appli/CMakeLists.txt new file mode 100644 index 0000000..71b459c --- /dev/null +++ b/appli/CMakeLists.txt @@ -0,0 +1,6 @@ + +SUBDIRS( + examples + ) + +## eof - $RCSfile$ diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt new file mode 100644 index 0000000..dff8fe0 --- /dev/null +++ b/appli/examples/CMakeLists.txt @@ -0,0 +1,31 @@ +IF(BUILD_EXAMPLES) + SET( + APPLIS + example_ImageAlgorithmRegionGrow_00 + example_ImageAlgorithmDijkstra_00 + example_ImageAlgorithmFastMarching_00 + ) + + FOREACH(APP ${APPLIS}) + ADD_EXECUTABLE(${APP} ${APP}.cxx) + TARGET_LINK_LIBRARIES(${APP} FrontAlgorithms) + ENDFOREACH(APP) + + IF(USE_VTK) + SET( + vtk_APPLIS + example_ImageAlgorithmRegionGrow_01 + example_ImageAlgorithmRegionGrow_MultipleThresholds + example_ImageAlgorithmDijkstra_01 + example_ImageAlgorithmDijkstra_02 + example_ImageAlgorithmFastMarching_01 + ) + + FOREACH(APP ${vtk_APPLIS}) + ADD_EXECUTABLE(${APP} ${APP}.cxx) + TARGET_LINK_LIBRARIES(${APP} FrontAlgorithms) + ENDFOREACH(APP) + ENDIF(USE_VTK) +ENDIF(BUILD_EXAMPLES) + +## eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmDijkstra_00.cxx b/appli/examples/example_ImageAlgorithmDijkstra_00.cxx new file mode 100644 index 0000000..86cff08 --- /dev/null +++ b/appli/examples/example_ImageAlgorithmDijkstra_00.cxx @@ -0,0 +1,44 @@ +#include +#include + +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef double TCost; +typedef itk::Image< TPixel, Dim > TImage; + +typedef fpa::Image::Dijkstra< TImage, TCost > TFrontAlgorithm; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Create a dummy image + TImage::SizeType imageSize; + imageSize.Fill( 100 ); + + TImage::SpacingType imageSpacing; + imageSpacing.Fill( 1 ); + + TImage::Pointer image = TImage::New( ); + image->SetRegions( imageSize ); + image->SetSpacing( imageSpacing ); + image->Allocate( ); + image->FillBuffer( TPixel( 1 ) ); + + // Seed + TImage::IndexType seed; + seed.Fill( 50 ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + algorithm->AddSeed( seed, TCost( 0 ) ); + algorithm->SetInput( image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->Update( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmDijkstra_01.cxx b/appli/examples/example_ImageAlgorithmDijkstra_01.cxx new file mode 100644 index 0000000..b2614ed --- /dev/null +++ b/appli/examples/example_ImageAlgorithmDijkstra_01.cxx @@ -0,0 +1,152 @@ +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + +typedef itk::ImageFileReader< TImage > TImageReader; +typedef fpa::Image::Dijkstra< TImage, TScalar > TFrontAlgorithm; + +typedef +fpa::VTK::Image2DObserver< TFrontAlgorithm, vtkRenderWindow > +TObserver; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image [stop_at_one_front]" << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + bool stop_at_one_front = false; + if( 2 < argc ) + stop_at_one_front = ( std::atoi( argv[ 2 ] ) == 1 ); + + // Read image + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + // VTK visualization + vtkSmartPointer< vtkImageActor > actor = + vtkSmartPointer< vtkImageActor >::New( ); + actor->SetInputData( vtk_image->GetOutput( ) ); + + vtkSmartPointer< vtkRenderer > renderer = + vtkSmartPointer< vtkRenderer >::New( ); + renderer->SetBackground( 0.1, 0.2, 0.7 ); + renderer->AddActor( actor ); + vtkSmartPointer< vtkRenderWindow > window = + vtkSmartPointer< vtkRenderWindow >::New( ); + window->SetSize( 800, 800 ); + window->AddRenderer( renderer ); + + // VTK interaction + vtkSmartPointer< vtkInteractorStyleImage > imageStyle = + vtkSmartPointer< vtkInteractorStyleImage >::New( ); + vtkSmartPointer< vtkRenderWindowInteractor > interactor = + vtkSmartPointer< vtkRenderWindowInteractor >::New( ); + interactor->SetInteractorStyle( imageStyle ); + window->SetInteractor( interactor ); + window->Render( ); + + // Create the widget and its representation + vtkSmartPointer< vtkPointHandleRepresentation3D > handle = + vtkSmartPointer< vtkPointHandleRepresentation3D >::New( ); + handle->GetProperty( )->SetColor( 1, 0, 0 ); + vtkSmartPointer< vtkSeedRepresentation > rep = + vtkSmartPointer< vtkSeedRepresentation >::New( ); + rep->SetHandleRepresentation( handle ); + + vtkSmartPointer< vtkSeedWidget > widget = + vtkSmartPointer< vtkSeedWidget >::New( ); + widget->SetInteractor( interactor ); + widget->SetRepresentation( rep ); + + // Let some interaction + interactor->Initialize( ); + window->Render( ); + widget->On( ); + interactor->Start( ); + + // Configure observer + TObserver::Pointer obs = TObserver::New( ); + obs->SetImage( input_image, window ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ ) + { + double pos[ 3 ]; + rep->GetSeedWorldPosition( s, pos ); + + TImage::PointType pnt; + pnt[ 0 ] = pos[ 0 ]; + pnt[ 1 ] = pos[ 1 ]; + + TImage::IndexType idx; + if( input_image->TransformPhysicalPointToIndex( pnt, idx ) ) + { + algorithm->AddSeed( idx, 0 ); + std::cout << " Seed --> " << idx << std::endl; + + } // fi + + } // rof + algorithm->AddObserver( itk::AnyEvent( ), obs ); + algorithm->ThrowEventsOn( ); + algorithm->SetInput( input_image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetStopAtOneFront( stop_at_one_front ); + algorithm->Update( ); + + // One last interaction + window->Render( ); + interactor->Start( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmDijkstra_02.cxx b/appli/examples/example_ImageAlgorithmDijkstra_02.cxx new file mode 100644 index 0000000..728313e --- /dev/null +++ b/appli/examples/example_ImageAlgorithmDijkstra_02.cxx @@ -0,0 +1,205 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef itk::RGBAPixel< TPixel > TRGBAPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TScalar, Dim > TDistanceMap; +typedef itk::Image< TRGBAPixel, Dim > TRGBAImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; +typedef itk::ImageToVTKImageFilter< TRGBAImage > TVTKRGBAImage; + +typedef itk::ImageFileReader< TImage > TImageReader; +typedef +itk::SignedDanielssonDistanceMapImageFilter< TImage, TDistanceMap > +TDistanceFilter; + +typedef fpa::Image::Dijkstra< TDistanceMap, TScalar > TFrontAlgorithm; +typedef fpa::Base::TreeExtractor< TFrontAlgorithm > TExtractor; +typedef fpa::Base::Functors::InvertCostFunction< TScalar > TScalarFunction; + +typedef +fpa::VTK::Image2DObserver< TExtractor, vtkRenderWindow > +TObserver; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image [stop_at_one_front]" << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + bool stop_at_one_front = false; + if( 2 < argc ) + stop_at_one_front = ( std::atoi( argv[ 2 ] ) == 1 ); + + // Read image + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + + TDistanceFilter::Pointer dist_filter = TDistanceFilter::New( ); + dist_filter->SetInput( input_image_reader->GetOutput( ) ); + dist_filter->InsideIsPositiveOn( ); + dist_filter->SquaredDistanceOff( ); + dist_filter->Update( ); + + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + // VTK visualization + vtkSmartPointer< vtkImageActor > actor = + vtkSmartPointer< vtkImageActor >::New( ); + actor->SetInputData( vtk_image->GetOutput( ) ); + + vtkSmartPointer< vtkRenderer > renderer = + vtkSmartPointer< vtkRenderer >::New( ); + renderer->SetBackground( 0.1, 0.2, 0.7 ); + renderer->AddActor( actor ); + vtkSmartPointer< vtkRenderWindow > window = + vtkSmartPointer< vtkRenderWindow >::New( ); + window->SetSize( 800, 800 ); + window->AddRenderer( renderer ); + + // VTK interaction + vtkSmartPointer< vtkInteractorStyleImage > imageStyle = + vtkSmartPointer< vtkInteractorStyleImage >::New( ); + vtkSmartPointer< vtkRenderWindowInteractor > interactor = + vtkSmartPointer< vtkRenderWindowInteractor >::New( ); + interactor->SetInteractorStyle( imageStyle ); + window->SetInteractor( interactor ); + window->Render( ); + + // Create the widget and its representation + vtkSmartPointer< vtkPointHandleRepresentation3D > handle = + vtkSmartPointer< vtkPointHandleRepresentation3D >::New( ); + handle->GetProperty( )->SetColor( 1, 0, 0 ); + vtkSmartPointer< vtkSeedRepresentation > rep = + vtkSmartPointer< vtkSeedRepresentation >::New( ); + rep->SetHandleRepresentation( handle ); + + vtkSmartPointer< vtkSeedWidget > widget = + vtkSmartPointer< vtkSeedWidget >::New( ); + widget->SetInteractor( interactor ); + widget->SetRepresentation( rep ); + + // Let some interaction + interactor->Initialize( ); + window->Render( ); + widget->On( ); + interactor->Start( ); + + // Configure observer + TObserver::Pointer obs = TObserver::New( ); + obs->SetImage( dist_filter->GetOutput( ), window ); + + // Configure membership function + TScalarFunction::Pointer cost_function = TScalarFunction::New( ); + + // Configure algorithm + TExtractor::Pointer algorithm = TExtractor::New( ); + for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ ) + { + double pos[ 3 ]; + rep->GetSeedWorldPosition( s, pos ); + + TImage::PointType pnt; + pnt[ 0 ] = pos[ 0 ]; + pnt[ 1 ] = pos[ 1 ]; + + TImage::IndexType idx; + if( input_image->TransformPhysicalPointToIndex( pnt, idx ) ) + { + algorithm->AddSeed( idx, 0 ); + if( s == 0 ) + algorithm->SetRoot( idx ); + else + algorithm->AddLeaf( idx ); + + } // fi + + } // rof + algorithm->AddObserver( itk::AnyEvent( ), obs ); + algorithm->ThrowEventsOn( ); + algorithm->SetInput( dist_filter->GetOutput( ) ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetStopAtOneFront( stop_at_one_front ); + algorithm->SetCostConversion( cost_function ); + algorithm->Update( ); + + // Create image with tree + TPixel transparent_color[] = { 0, 0, 0, 0 }; + TPixel solid_color[] = { 255, 255, 0, 255 }; + TRGBAImage::Pointer tree_image = TRGBAImage::New( ); + tree_image->CopyInformation( input_image ); + tree_image->SetRequestedRegionToLargestPossibleRegion( ); + tree_image->SetBufferedRegion( tree_image->GetRequestedRegion( ) ); + tree_image->Allocate( ); + tree_image->FillBuffer( TRGBAPixel( transparent_color ) ); + + const TExtractor::TTree& tree = algorithm->GetTree( ); + TExtractor::TTree::const_iterator tIt = tree.begin( ); + for( ; tIt != tree.end( ); ++tIt ) + tree_image->SetPixel( tIt->first, TRGBAPixel( solid_color ) ); + + TVTKRGBAImage::Pointer vtk_tree_image = TVTKRGBAImage::New( ); + vtk_tree_image->SetInput( tree_image ); + vtk_tree_image->Update( ); + + vtkSmartPointer< vtkImageActor > tree_actor = + vtkSmartPointer< vtkImageActor >::New( ); + tree_actor->SetInputData( vtk_tree_image->GetOutput( ) ); + + renderer->AddActor( tree_actor ); + window->Render( ); + + // One last interaction + interactor->Start( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmFastMarching_00.cxx b/appli/examples/example_ImageAlgorithmFastMarching_00.cxx new file mode 100644 index 0000000..1573ce6 --- /dev/null +++ b/appli/examples/example_ImageAlgorithmFastMarching_00.cxx @@ -0,0 +1,44 @@ +#include +#include + +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef double TCost; +typedef itk::Image< TPixel, Dim > TImage; + +typedef fpa::Image::FastMarching< TImage, TCost > TFrontAlgorithm; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Create a dummy image + TImage::SizeType imageSize; + imageSize.Fill( 100 ); + + TImage::SpacingType imageSpacing; + imageSpacing.Fill( 1 ); + + TImage::Pointer image = TImage::New( ); + image->SetRegions( imageSize ); + image->SetSpacing( imageSpacing ); + image->Allocate( ); + image->FillBuffer( TPixel( 1 ) ); + + // Seed + TImage::IndexType seed; + seed.Fill( 50 ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + algorithm->AddSeed( seed, TCost( 0 ) ); + algorithm->SetInput( image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->Update( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmFastMarching_01.cxx b/appli/examples/example_ImageAlgorithmFastMarching_01.cxx new file mode 100644 index 0000000..cfa306f --- /dev/null +++ b/appli/examples/example_ImageAlgorithmFastMarching_01.cxx @@ -0,0 +1,152 @@ +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + +typedef itk::ImageFileReader< TImage > TImageReader; +typedef fpa::Image::FastMarching< TImage, TScalar > TFrontAlgorithm; + +typedef +fpa::VTK::Image2DObserver< TFrontAlgorithm, vtkRenderWindow > +TObserver; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image [stop_at_one_front]" << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + bool stop_at_one_front = false; + if( 2 < argc ) + stop_at_one_front = ( std::atoi( argv[ 2 ] ) == 1 ); + + // Read image + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + // VTK visualization + vtkSmartPointer< vtkImageActor > actor = + vtkSmartPointer< vtkImageActor >::New( ); + actor->SetInputData( vtk_image->GetOutput( ) ); + + vtkSmartPointer< vtkRenderer > renderer = + vtkSmartPointer< vtkRenderer >::New( ); + renderer->SetBackground( 0.1, 0.2, 0.7 ); + renderer->AddActor( actor ); + vtkSmartPointer< vtkRenderWindow > window = + vtkSmartPointer< vtkRenderWindow >::New( ); + window->SetSize( 800, 800 ); + window->AddRenderer( renderer ); + + // VTK interaction + vtkSmartPointer< vtkInteractorStyleImage > imageStyle = + vtkSmartPointer< vtkInteractorStyleImage >::New( ); + vtkSmartPointer< vtkRenderWindowInteractor > interactor = + vtkSmartPointer< vtkRenderWindowInteractor >::New( ); + interactor->SetInteractorStyle( imageStyle ); + window->SetInteractor( interactor ); + window->Render( ); + + // Create the widget and its representation + vtkSmartPointer< vtkPointHandleRepresentation3D > handle = + vtkSmartPointer< vtkPointHandleRepresentation3D >::New( ); + handle->GetProperty( )->SetColor( 1, 0, 0 ); + vtkSmartPointer< vtkSeedRepresentation > rep = + vtkSmartPointer< vtkSeedRepresentation >::New( ); + rep->SetHandleRepresentation( handle ); + + vtkSmartPointer< vtkSeedWidget > widget = + vtkSmartPointer< vtkSeedWidget >::New( ); + widget->SetInteractor( interactor ); + widget->SetRepresentation( rep ); + + // Let some interaction + interactor->Initialize( ); + window->Render( ); + widget->On( ); + interactor->Start( ); + + // Configure observer + TObserver::Pointer obs = TObserver::New( ); + obs->SetImage( input_image, window ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ ) + { + double pos[ 3 ]; + rep->GetSeedWorldPosition( s, pos ); + + TImage::PointType pnt; + pnt[ 0 ] = pos[ 0 ]; + pnt[ 1 ] = pos[ 1 ]; + + TImage::IndexType idx; + if( input_image->TransformPhysicalPointToIndex( pnt, idx ) ) + { + algorithm->AddSeed( idx, 0 ); + std::cout << " Seed --> " << idx << std::endl; + + } // fi + + } // rof + algorithm->AddObserver( itk::AnyEvent( ), obs ); + algorithm->ThrowEventsOn( ); + algorithm->SetInput( input_image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetStopAtOneFront( stop_at_one_front ); + algorithm->Update( ); + + // One last interaction + window->Render( ); + interactor->Start( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmRegionGrow_00.cxx b/appli/examples/example_ImageAlgorithmRegionGrow_00.cxx new file mode 100644 index 0000000..149977a --- /dev/null +++ b/appli/examples/example_ImageAlgorithmRegionGrow_00.cxx @@ -0,0 +1,51 @@ +#include +#include + +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef itk::Image< TPixel, Dim > TImage; + +typedef fpa::Image::RegionGrow< TImage > TFrontAlgorithm; +typedef +fpa::Image::Functors::RegionGrowAllBelongsFunction< TImage > +TMembershipFunction; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Create a dummy image + TImage::SizeType imageSize; + imageSize.Fill( 10 ); + + TImage::SpacingType imageSpacing; + imageSpacing.Fill( 1 ); + + TImage::Pointer image = TImage::New( ); + image->SetRegions( imageSize ); + image->SetSpacing( imageSpacing ); + image->Allocate( ); + image->FillBuffer( TPixel( 1 ) ); + + // Seed + TImage::IndexType seed; + seed.Fill( 5 ); + + // Configure membership function + TMembershipFunction::Pointer membership = TMembershipFunction::New( ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + algorithm->AddSeed( seed, 0 ); + algorithm->SetInput( image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetMembershipFunction( membership ); + algorithm->Update( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmRegionGrow_01.cxx b/appli/examples/example_ImageAlgorithmRegionGrow_01.cxx new file mode 100644 index 0000000..e41b65d --- /dev/null +++ b/appli/examples/example_ImageAlgorithmRegionGrow_01.cxx @@ -0,0 +1,157 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 2; +typedef unsigned char TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + +typedef itk::ImageFileReader< TImage > TImageReader; +typedef fpa::Image::RegionGrow< TImage > TFrontAlgorithm; + +typedef +fpa::Image::Functors::RegionGrowThresholdFunction< TImage > +TMembershipFunction; + +typedef +fpa::VTK::Image2DObserver< TFrontAlgorithm, vtkRenderWindow > +TObserver; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 4 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image lower_threshold upper_threshold" << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + TPixel lower_threshold = TPixel( std::atof( argv[ 2 ] ) ); + TPixel upper_threshold = TPixel( std::atof( argv[ 3 ] ) ); + + // Read image + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + // VTK visualization + vtkSmartPointer< vtkImageActor > actor = + vtkSmartPointer< vtkImageActor >::New( ); + actor->SetInputData( vtk_image->GetOutput( ) ); + + vtkSmartPointer< vtkRenderer > renderer = + vtkSmartPointer< vtkRenderer >::New( ); + renderer->SetBackground( 0.1, 0.2, 0.7 ); + renderer->AddActor( actor ); + vtkSmartPointer< vtkRenderWindow > window = + vtkSmartPointer< vtkRenderWindow >::New( ); + window->SetSize( 800, 800 ); + window->AddRenderer( renderer ); + + // VTK interaction + vtkSmartPointer< vtkInteractorStyleImage > imageStyle = + vtkSmartPointer< vtkInteractorStyleImage >::New( ); + vtkSmartPointer< vtkRenderWindowInteractor > interactor = + vtkSmartPointer< vtkRenderWindowInteractor >::New( ); + interactor->SetInteractorStyle( imageStyle ); + window->SetInteractor( interactor ); + window->Render( ); + + // Create the widget and its representation + vtkSmartPointer< vtkPointHandleRepresentation3D > handle = + vtkSmartPointer< vtkPointHandleRepresentation3D >::New( ); + handle->GetProperty( )->SetColor( 1, 0, 0 ); + vtkSmartPointer< vtkSeedRepresentation > rep = + vtkSmartPointer< vtkSeedRepresentation >::New( ); + rep->SetHandleRepresentation( handle ); + + vtkSmartPointer< vtkSeedWidget > widget = + vtkSmartPointer< vtkSeedWidget >::New( ); + widget->SetInteractor( interactor ); + widget->SetRepresentation( rep ); + + // Let some interaction + interactor->Initialize( ); + window->Render( ); + widget->On( ); + interactor->Start( ); + + // Configure observer + TObserver::Pointer obs = TObserver::New( ); + obs->SetImage( input_image, window ); + + // Configure membership function + TMembershipFunction::Pointer membership = TMembershipFunction::New( ); + membership->SetLowerThreshold( lower_threshold ); + membership->SetUpperThreshold( upper_threshold ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ ) + { + double pos[ 3 ]; + rep->GetSeedWorldPosition( s, pos ); + + TImage::PointType pnt; + pnt[ 0 ] = TScalar( pos[ 0 ] ); + pnt[ 1 ] = TScalar( pos[ 1 ] ); + + TImage::IndexType idx; + if( input_image->TransformPhysicalPointToIndex( pnt, idx ) ) + algorithm->AddSeed( idx, 0 ); + + } // rof + algorithm->AddObserver( itk::AnyEvent( ), obs ); + algorithm->ThrowEventsOn( ); + algorithm->SetInput( input_image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetMembershipFunction( membership ); + algorithm->Update( ); + + // One last interaction + interactor->Start( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx new file mode 100644 index 0000000..165dafb --- /dev/null +++ b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx @@ -0,0 +1,118 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef short TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + +typedef itk::ImageFileReader< TImage > TImageReader; + +typedef +fpa::Image::RegionGrowWithMultipleThresholds< TImage > +TFrontAlgorithm; +typedef +fpa::VTK::Image3DObserver< TFrontAlgorithm, vtkRenderWindow > +TObserver; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 8 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image s_x s_y s_z thr_0 thr_1 n_samples" << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + TImage::PointType seed_pnt; + seed_pnt[ 0 ] = std::atof( argv[ 2 ] ); + seed_pnt[ 1 ] = std::atof( argv[ 3 ] ); + seed_pnt[ 2 ] = std::atof( argv[ 4 ] ); + TPixel thr_0 = TPixel( std::atof( argv[ 5 ] ) ); + TPixel thr_1 = TPixel( std::atof( argv[ 6 ] ) ); + unsigned int n_samples = std::atoi( argv[ 7 ] ); + + // Read image + TImageReader::Pointer input_image_reader = TImageReader::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + TImage::SpacingType spac = input_image->GetSpacing( ); + double min_spac = spac[ 0 ]; + for( unsigned int d = 1; d < Dim; d++ ) + min_spac = ( spac[ d ] < min_spac )? spac[ d ]: min_spac; + + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + vtkSmartPointer< vtkSphereSource > seed = + vtkSmartPointer< vtkSphereSource >::New( ); + seed->SetCenter( seed_pnt[ 0 ], seed_pnt[ 1 ], seed_pnt[ 2 ] ); + seed->SetRadius( min_spac * double( 5 ) ); + seed->Update( ); + + fpa::VTK::ImageMPR view; + view.SetBackground( 0.3, 0.2, 0.8 ); + view.SetSize( 800, 800 ); + view.SetImage( vtk_image->GetOutput( ) ); + view.AddPolyData( seed->GetOutput( ), 1, 0, 0 ); + view.Start( ); + + // Configure observer + TObserver::Pointer obs = TObserver::New( ); + obs->SetImage( input_image, view.GetWindow( ) ); + + // Configure algorithm + TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( ); + algorithm->AddThresholds( thr_0, thr_1, n_samples ); + + TImage::IndexType seed_idx; + input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); + algorithm->AddSeed( seed_idx, 0 ); + + algorithm->AddObserver( itk::AnyEvent( ), obs ); + algorithm->ThrowEventsOn( ); + algorithm->SetInput( input_image ); + algorithm->SetNeighborhoodOrder( 1 ); + algorithm->SetDerivativeThreshold( double( 3 ) ); + algorithm->Update( ); + + view.Start( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt new file mode 100644 index 0000000..0123559 --- /dev/null +++ b/cmake/CMakeLists.txt @@ -0,0 +1,7 @@ +CONFIGURE_FILE( + FrontAlgorithmsConfig.cmake.in + ${PROJECT_BINARY_DIR}/FrontAlgorithmsConfig.cmake + @ONLY + ) + +## eof - $RCSfile$ diff --git a/cmake/FrontAlgorithmsConfig.cmake.in b/cmake/FrontAlgorithmsConfig.cmake.in new file mode 100644 index 0000000..da3d3e8 --- /dev/null +++ b/cmake/FrontAlgorithmsConfig.cmake.in @@ -0,0 +1,43 @@ +FIND_PATH( + FrontAlgorithms_INCLUDE_DIR1 + fpa/Base/Algorithm.h + PATHS + /usr/include + /usr/local/include + @PROJECT_SOURCE_DIR@/lib + @PROJECT_BINARY_DIR@/lib + @CMAKE_INSTALL_PREFIX@/include + ) + +FIND_PATH( + FrontAlgorithms_INCLUDE_DIR2 + fpa/FrontAlgorithms_Export.h + PATHS + /usr/include + /usr/local/include + @PROJECT_SOURCE_DIR@/lib + @PROJECT_BINARY_DIR@/lib + @CMAKE_INSTALL_PREFIX@/include + ) + +INCLUDE_DIRECTORIES( + ${FrontAlgorithms_INCLUDE_DIR1} + ${FrontAlgorithms_INCLUDE_DIR2} + ) + +FIND_LIBRARY( + FrontAlgorithms_LIBRARY_NAME + FrontAlgorithms + PATHS + /usr/lib + /usr/local/lib + @PROJECT_BINARY_DIR@ + @CMAKE_INSTALL_PREFIX@/lib + ) + +SET( + FrontAlgorithms_LIBRARIES + ${FrontAlgorithms_LIBRARY_NAME} + ) + +## eof - $RCSfile$ diff --git a/data/binary_test_2D_00.png b/data/binary_test_2D_00.png new file mode 100644 index 0000000000000000000000000000000000000000..fbe9d39bd59189ecf7f8ef87e48abde2f3d01953 GIT binary patch literal 5192 zcma)9c{o)6_aD=ZI}EO|PKi5=td%v|KIR(BAZe2=Yhx|DN_LmVP#8-{ib^R-krY`A z@reqNiYy^ZWhrGDe6K#w@BiO@p1Yj$exLI?=XK6=-_J=nc)&!Ee-l3rhZ8hq7#_yq zHo)sQh$(QZn{CUn7r~!l8-&9Nh_2r}xSU)mOvx8)YGK6pnO|H;NsAJfXoRV@cpbF~ zru+H$xcLTSA`V9nbaM%IW5Q>>f<0jqQ;UPl?QVQHoM^47A#yb2`nC({b$L%f zAN+ba&2J>U#}W!JiLbGyhW19FGcQ&}0n5ef8#3Z$LummE*9jo%!F{{Y?1ivTz9#oA zAJ8Bg=?R=3ia&5rBGU7GSY-yz{}c%I=To_Fh51pwmA4}9kaW2vnja-FMQ9p~zkN_3 za-9FxV0`UCuv*_>is&{N-?N}t-F5-5X$EP3vu(f53v|e+)C&u`^2V6P!Om>J=>9Q> zet8n0xh^fdD_k-}X>2|6!AB;sFwACAK&1R)H_5QeMpZw!Jv~HRpvLXAM84}0UoU9_ zo-z#gen_45)Z#>zHu~B_pt2C^JJRHBx0UE>r-x8gGOoEGicVJbK`=pHy?g-8X4 zdR5-9wspO|9AQBY40fx$S8MyDXX!#bgC%WmQ&Qx2ENL{v{ffjUk^*w^mmEkJvO*~qK*vD;16xxl$&Ey=i`C^Hj<_BqHk;-bt13Yk zx^^TWmTC&I_SPoG8gb-x1&04H88-u#t(UFTSnA9)xPe2F2kD_yMiW2Z&`xrom7fgx zIR_y;hgfYN>!#u<&MjeDAbre&aYe8)Qfh7n|7fjKa&zEHusT4zEY|v>Zcx=Je4+U8 zol`Y2c~Kpcc^<{j>{G?e3oScql{sr`bx{Kcgje$4F8B47B_1!+j9QMWm5H$Rh^lfd zEjxU9cVtZOd#CbE3id@66#I={o{Y#9+H|$#RPyGmJl5dnYpuIcDKk4Fe^)c7dMM#@ znC|rZyrH%&4mP`V;v+<6{U*b7iu~T2_vF8F(Jp$i?gMIBfkOL1Deo29SrMiQYs zcHMK5(E0nNh5N5tCVe{o!iP?87}zuDyL_Yjo>h!dOLF$dT$pwXw^*njMuSAkbPbk# z2IjlkI2KtKbFF)+E!n&L0<422XPbA}>BQD6+{=YO$;DB8D#fH8_&Spe89SsLDx?Fl z@L&1`_oJW}ivewPZ$UTw@U|KJDV}^taEh&}!WJ$%tM7Rta+I47J=okik$n2J+CvQo zmYPUo!RML`*_C@(mDHq=ozA>5pz=00$b(498H5_Q%)fxNs{Xx9rY6R-L<(r^#vlR2RCs)+Xj9uq5g=V|xDWI$)ta_3e3 z{Pq$whojqnOsQU>0S$jbnNdvy{Lhzs)?sKqOV&t79!%8!EkZ|O|GG5zqi=nzj7=~E z7-qs!LocWFEeL4ttFeuILECL~k-u2vG!g88q}i*Xu583-z^wXr|Qd)5wOhdrmYgqK09!vGCx_@481& zDI{?zRY8UWCr>I%hywfWC&{m=AaA$oFByZ#%(p`_0b>I-%I`o_w>hOk#DHTLJboC6q?an1%Mc)Jo`~dnk)hhItr)3HbLnVk`$n!zwpnUfp0|9cIg>FxWyv`giiTGS`r15iCd$dr?4P}H0mc_h zYUCzeTlwHo4VUfgkh+nmy_H$#Nu{lD{9HlB1}ZW35Cy{0s#6D40i_-4=e zlaBv(RLiCn5i}KO<%8DCq%Hi)$iqoy0xzm2%sFSFcNRZ!2`;T9Y>_zW z&N0oyciaYeWOggmE+29c5^sD777nZ$w}vMu&4EBzKWVq%*pc6n?0rTh9KG!fEfG6ObrYsRh7e^Ys$mt^S(+7FLyk0 z9DR}>BO8fis%CPf^871L51j^P&sRhp5qRH{R@kc(?65~@=Dx$S>*JEiZG+CIU+lhn zku8{Is5^tHG?jag?~QwsNsKb_ega3N@V&A#gHHA`%!0SzpqD;+m>+6Kv_03 z6Q(&J)p1sfkzGD#=W0AI!M;1;mUpRE)xsi;5^M-{KIC{cJX+F(^Do-Sx^^F;cwJEd zCHNQ-$eD1Pc3lafwmk%Nw$OUN7`3822~>7aCzQ<%Kb-khq{i_5fEO`?%^^8^teB%-w(WZewqj*uAWX;LxX< zN>WI<*xDyT+DA#C!wF{Nt}oQqsK-Wg?*jyfUt> z9M7+6zVErZyPB~T8DtqjL& zxvpx!o3Vyqil615zH;5B-s&!ekD7`Ulce&^AF34mY1aHK>I%T~2^>y) z>@K3KKehpV>Oy`eIH7lf;m_@VS$%sukyLRTjs%BsQb)(k*KDrK)hV%YUM`H}U&}G^ z)!M?hwd+M4My!f271D~PDM(q^R7-0~%4(au8Pw)E9WPcpRy8gCIrv3rQ8}rgdnc>2 z`gr!qhg?md8<;%;v2yS14;C?e)pp#sIAX6B+x?_dwxVk2grw@RH_y?yEk;A$j-gFO z_F1DtcAr(*lBMpONMW+y4=)-E2wR(N20$vZZKDVuVkg`e|2O9bkF2N8yX$JA&Q4oRv*K`;U6q=bs z9?GB`P3%uOd|)>Ndg}7sP8aEWn2`kw9fkm`MR}0usHmM(q(U$iZ}>MFy!=viw`jo( zi2hXA%9a-cVm6@g(gq+#9))=sJaYd9eOrhEbZmwhb?XeQvtkUP^`HT?o(~>syqK(L zvJWFeJO>LqnhL^5+xQ(33ZQW_c3tu5%Ti7mB+oWuU`$}L#M!fGJ!Ak`@SxvX!m+@> zw}Z9k+HLvORgj9~?MX1e!}?sh4s$_{O3(}iW?kbk;&Y{dyAm)5re)D54j{*`ZDa!- zd(-{;{F?%KqhY8B`|%0;qGF$v%)$oLzf@G)~)9ugocQ(DTsOu&+ml_@J-xoj!m z3M%*el&)S@0O?f=0V9_Q&pA6Ib3@!w^@0>xzH5u24guN|6*Fz5ILi0rs;s=llD!e9P8zm>g4%b%D`7{5n~`8`wr;k zhyJRD|8X_qU|Ea^t+RQjVS^u~onN?Pu4H7Tn)El3Xn4A>z@w&qV$=?^h`}|8q9Kn- zK)@vCfei(U7@o|q1?lV7a&hnzAPaxpn}7IQqYRQ+lYB{GO{J~M)k_dEjMh^G z=e~HHhuM32wBBqpI0lkV+R^EvfZtNXQ*ks7a|N)?HzS51i0+QB2|}@)a8fmalQ7@{`WtYZ+5C{9a4eiNJ!(=gz>B#0 zMtM4;(Xtm3hLglTt;{@2r!wAa*+oJlT#amQJr^K{92$%#H%c9@7pAkx0Nhgi?4p1K z(kB)zWLRH z+ zrL!So8DGw(T#7iG1e+6E)6ZZ1-ChUW4K+`6J;VK-AZ5VG+Vr^hy>H`}lzNU4n!8YQ z-u8OYN}Ui2YYDc>8HrQ6E63DE|FmaXarGi@=eO_5*s*=9z@>Kz7vOH9vbgNx3kc27 zEa;t2qyWlpjLO+q3UK!XL!x*q9X#B~QL}C1x`2p^RGQJrRtP<+PQ7jsC$0{MytwzhR%xE5SrbX_&b%sc3QKH}YoJKA);Y_j&TNSiM&a5gm6dXBE zIr+QFKhAJSL5rCLTL~(!&UX=QI*2Psp7vf^Zg`Uls zQzwZ}9^3s;iNgr6rH`$hkz}1xwJ%Yjlj7i)MBCZca8~g>yAl;V9Sk&1)yXSjw$y#~ z{&7}qZUbp#CFpdwVsq}sJ?5GAnSPjuD_0-Q)T|Oe=){IUo-SoxJveLmy65V2S3&Q` z(9z`c*c9^S&+9Q=NSzt|Xecd?0eJ?UH%cXT;7NdgI2+mrV`n1`J2t*gy8|&G27^cIw_JebuptMN3_%2E zBDAS5^H&TIFdPFrgZ}~jUtE7rkrV%&_v>bZqDcNI0k7xjpZjhj3Bbl1Jp=8S^(pV4 zL}*{W9)=4vD6+*qXUK|2F~dJ0LZbRVE_OC literal 0 HcmV?d00001 diff --git a/data/ones_image.png b/data/ones_image.png new file mode 100644 index 0000000000000000000000000000000000000000..f7afaab3a8241507fc085fd4450f6bfa393ef4db GIT binary patch literal 1697 zcmeAS@N?(olHy`uVBq!ia0y~yU{(NO4xj+TKk3LVK#H@#BeIx*fm;}a85w5HkpK#^ zmw5WRvftxo=F>J^eIwQZC?uO15>euupPQSSR|4cRFgO>bCYGe8D3oWGWGJ|M`Ua%v zrLr?Hu!(!RIEGZ*dV6^zF9QR|VS`iu-M=%HvHfoG@MPPuh9yD3fJ1_Xhl!1`xxw*( zLV^GY6B2@`L{Wty!~iu9Nfm|=1Kd`)Dl9@Cm|{roh3Ui)Lh%(u6^0N@Gm={|gfN1` zfCE#gV>BR<6A>s7M*|X(sUTs2k^Fd=FoZBdZ8RVmN3%A#WQHdV%<{0A#4_?jG~?Ve VY>s=^oN@!z9iFa!F6*2UngGDBTyX#Z literal 0 HcmV?d00001 diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt new file mode 100644 index 0000000..f1821ee --- /dev/null +++ b/lib/CMakeLists.txt @@ -0,0 +1,67 @@ +SET(LIB_NAME FrontAlgorithms) + +## ================ +## = Source files = +## ================ + +CONFIGURE_FILE( + fpa/Common.cxx.in + ${PROJECT_BINARY_DIR}/lib/fpa/Common.cxx + ) + +FILE(GLOB ${LIB_NAME}_HEADERS "fpa/*.h" "fpa/*.hxx") +FILE(GLOB ${LIB_NAME}_BASE_HEADERS "fpa/Base/*.h" "fpa/Base/*.hxx") +FILE(GLOB ${LIB_NAME}_IMAGE_HEADERS "fpa/Image/*.h" "fpa/Image/*.hxx") + +FILE(GLOB ${LIB_NAME}_SOURCES "fpa/*.cxx") +FILE(GLOB ${LIB_NAME}_BASE_SOURCES "fpa/Base/*.cxx") +FILE(GLOB ${LIB_NAME}_IMAGE_SOURCES "fpa/Image/*.cxx") + +IF(USE_VTK) + FILE(GLOB ${LIB_NAME}_VTK_HEADERS "fpa/VTK/*.h" "fpa/VTK/*.hxx") + FILE(GLOB ${LIB_NAME}_VTK_SOURCES "fpa/VTK/*.cxx") +ENDIF(USE_VTK) + +SET( + ${LIB_NAME}_ALL_SOURCES + ${PROJECT_BINARY_DIR}/lib/fpa/Common.cxx + ${${LIB_NAME}_SOURCES} + ${${LIB_NAME}_BASE_SOURCES} + ${${LIB_NAME}_IMAGE_SOURCES} + ${${LIB_NAME}_VTK_SOURCES} + ) + +## ============= +## = Libraries = +## ============= + +SET( + ${LIB_NAME}_LINK_LIBRARIES + ${ITK_LIBRARIES} + ${VTK_LIBRARIES} + ) + +## ===================== +## = Compilation rules = +## ===================== + +ADD_LIBRARY( + ${LIB_NAME} + ${LIB_TYPE} + ${${LIB_NAME}_ALL_SOURCES} + ) +GENERATE_EXPORT_HEADER( + ${LIB_NAME} + BASE_NAME ${LIB_NAME} + EXPORT_MACRO_NAME ${LIB_NAME}_EXPORT + EXPORT_FILE_NAME ${PROJECT_BINARY_DIR}/lib/fpa/${LIB_NAME}_Export.h + STATIC_DEFINE ${LIB_NAME}_BUILT_AS_STATIC + ) +TARGET_LINK_LIBRARIES( + ${LIB_NAME} + ${${LIB_NAME}_LINK_LIBRARIES} + ${ITK_LIBRARIES} + ${VTK_LIBRARIES} + ) + +## eof - $RCSfile$ diff --git a/lib/fpa/Base/Algorithm.h b/lib/fpa/Base/Algorithm.h new file mode 100644 index 0000000..187a281 --- /dev/null +++ b/lib/fpa/Base/Algorithm.h @@ -0,0 +1,153 @@ +#ifndef __FPA__BASE__ALGORITHM__H__ +#define __FPA__BASE__ALGORITHM__H__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + * Base front propagation algorithm. From a series of start seeds with + * costs, a priority queue is filled and emptied updating costs. Each + * vertex could be marked as "visited", "in the front", "not yet there" + * or "freezed". + * + * @param T Traits used for this algorithm + */ + template< class T, class B > + class Algorithm + : public B + { + public: + // Standard class typdedefs + typedef Algorithm Self; + typedef B Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /// Templated types + typedef T TTraits; + typedef B TBaseFilter; + typedef typename T::TCost TCost; + typedef typename T::TResult TResult; + typedef typename T::TVertex TVertex; + typedef typename T::TVertexValue TVertexValue; + + protected: + typedef typename T::TFrontId _TFrontId; + typedef typename T::TNode _TNode; + typedef typename T::TNodes _TNodes; + typedef typename T::TVertexCmp _TVertexCmp; + + typedef std::map< TVertex, _TNode, _TVertexCmp > _TMarks; + + typedef std::pair< TVertex, bool > _TCollision; + typedef std::vector< _TCollision > _TCollisionSitesRow; + typedef std::vector< _TCollisionSitesRow > _TCollisionSites; + + public: + typedef BaseEvent< _TNode > TEvent; + typedef FrontEvent< _TNode > TFrontEvent; + typedef MarkEvent< _TNode > TMarkEvent; + typedef CollisionEvent< _TNode > TCollisionEvent; + typedef EndEvent< _TNode > TEndEvent; + + public: + itkTypeMacro( Algorithm, itkProcessObject ); + + itkBooleanMacro( StopAtOneFront ); + itkBooleanMacro( ThrowEvents ); + + itkGetConstMacro( StopAtOneFront, bool ); + itkGetConstMacro( ThrowEvents, bool ); + + itkSetMacro( StopAtOneFront, bool ); + itkSetMacro( ThrowEvents, bool ); + + public: + virtual void InvokeEvent( const itk::EventObject& e ); + virtual void InvokeEvent( const itk::EventObject& e ) const; + + /// Seeds manipulation + void AddSeed( const TVertex& s, const TResult& v ); + const TVertex& GetSeed( const unsigned long& i ) const; + void ClearSeeds( ); + unsigned long GetNumberOfSeeds( ) const; + + protected: + Algorithm( ); + virtual ~Algorithm( ); + + /// itk::ProcessObject + virtual void GenerateData( ); + + /// Base interface + virtual void _BeforeLoop ( ); + virtual void _AfterLoop ( ); + virtual void _Loop ( ); + virtual bool _CheckCollisions ( const _TNode& a, const _TNode& b ); + virtual bool _CheckStopCondition ( ); + virtual bool _UpdateResult ( _TNode& n ); + + /// Marks management + virtual void _InitializeMarks ( ); + virtual bool _IsMarked ( const TVertex& v ) const; + virtual _TFrontId _FrontId ( const TVertex& v ) const; + virtual TVertex _Parent ( const TVertex& v ) const; + virtual void _Mark ( const _TNode& n ); + + /// Pure virtual interface: vertices + virtual unsigned long _NumberOfVertices ( ) const = 0; + virtual TVertexValue _Value ( const TVertex& v ) const = 0; + virtual TResult _Result ( const TVertex& v ) const = 0; + + /// Pure virtual interface: edges + virtual double _Norm ( const TVertex& a, const TVertex& b ) const = 0; + virtual bool _Edge ( const TVertex& a, const TVertex& b ) const = 0; + virtual TCost _Cost ( const TVertex& a, const TVertex& b ) const = 0; + + /// Pure virtual interface: neighborhood + virtual void _Neighs ( const _TNode& n, _TNodes& N ) const = 0; + virtual bool _UpdateNeigh ( _TNode& nn, const _TNode& n ) = 0; + virtual void _NeighsInDim ( const _TNode& n, + const unsigned int& d, + _TNodes& N ) = 0; + + /// Pure virtual interface: queue + virtual void _InitializeQueue ( ) = 0; + virtual bool _IsQueueEmpty ( ) const = 0; + virtual void _QueuePush ( const _TNode& n ) = 0; + virtual _TNode _QueuePop ( ) = 0; + virtual void _QueueClear ( ) = 0; + + /// Pure virtual interface: results + virtual void _InitializeResults ( ) = 0; + + private: + // Purposely not implemented + Algorithm( const Self& ); + void operator=( const Self& ); + + protected: + bool m_StopAtOneFront; + bool m_ThrowEvents; + + _TNodes m_Seeds; + _TMarks m_Marks; + + _TCollisionSites m_CollisionSites; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__BASE__ALGORITHM__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Algorithm.hxx b/lib/fpa/Base/Algorithm.hxx new file mode 100644 index 0000000..286343f --- /dev/null +++ b/lib/fpa/Base/Algorithm.hxx @@ -0,0 +1,300 @@ +#ifndef __FPA__BASE__ALGORITHM__HXX__ +#define __FPA__BASE__ALGORITHM__HXX__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +InvokeEvent( const itk::EventObject& e ) +{ + if( this->m_ThrowEvents ) + this->Superclass::InvokeEvent( e ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +InvokeEvent( const itk::EventObject& e ) const +{ + if( this->m_ThrowEvents ) + this->Superclass::InvokeEvent( e ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +AddSeed( const TVertex& s, const TResult& v ) +{ + this->m_Seeds.push_back( _TNode( s, v, this->m_Seeds.size( ) ) ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +const typename fpa::Base::Algorithm< T, B >:: +TVertex& fpa::Base::Algorithm< T, B >:: +GetSeed( const unsigned long& i ) const +{ + return( this->m_Seeds[ i ].Vertex ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +ClearSeeds( ) +{ + this->m_Seeds.clear( ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +unsigned long fpa::Base::Algorithm< T, B >:: +GetNumberOfSeeds( ) const +{ + return( this->m_Seeds.size( ) ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +fpa::Base::Algorithm< T, B >:: +Algorithm( ) + : Superclass( ), + m_StopAtOneFront( false ), + m_ThrowEvents( false ) +{ +} + +// ------------------------------------------------------------------------- +template< class T, class B > +fpa::Base::Algorithm< T, B >:: +~Algorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +GenerateData( ) +{ + this->AllocateOutputs( ); + + unsigned long N = this->m_Seeds.size( ); + if( N == 0 ) + return; + this->m_CollisionSites.clear( ); + this->m_CollisionSites. + resize( N, _TCollisionSitesRow( N, _TCollision( TVertex( ), false ) ) ); + this->_BeforeLoop( ); + this->_Loop( ); + this->_AfterLoop( ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +_BeforeLoop( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +_AfterLoop( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +_Loop( ) +{ + this->_InitializeMarks( ); + this->_InitializeResults( ); + this->_InitializeQueue( ); + while( !( this->_IsQueueEmpty( ) ) ) + { + _TNode n = this->_QueuePop( ); + if( this->_IsMarked( n.Vertex ) ) + continue; + + this->_Mark( n ); + this->InvokeEvent( TMarkEvent( n ) ); + if( this->_UpdateResult( n ) ) + { + if( !( this->_CheckStopCondition( ) ) ) + { + _TNodes N; + this->_Neighs( n, N ); + typename _TNodes::iterator nnIt = N.begin( ); + while( nnIt != N.end( ) ) + { + if( this->_IsMarked( nnIt->Vertex ) ) + { + // Update real front identifier + nnIt->FrontId = this->_FrontId( nnIt->Vertex ); + + // Update collisions + if( this->_CheckCollisions( n, *nnIt ) ) + { + if( this->m_StopAtOneFront ) + { + this->_QueueClear( ); + nnIt = N.end( ); + + } // fi + + } // fi + } + else + { + if( this->_UpdateNeigh( *nnIt, n ) ) + { + nnIt->Parent = n.Vertex; + this->_QueuePush( *nnIt ); + this->InvokeEvent( TFrontEvent( *nnIt ) ); + + } // fi + + } // fi + if( nnIt != N.end( ) ) + nnIt++; + + } // elihw + } + else + this->_QueueClear( ); + + } // fi + + } // elihw + this->InvokeEvent( TEndEvent( ) ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +bool fpa::Base::Algorithm< T, B >:: +_CheckCollisions( const _TNode& a, const _TNode& b ) +{ + bool ret = false; + if( a.FrontId != b.FrontId ) + { + // Mark collision, if it is new + bool exists = this->m_CollisionSites[ a.FrontId ][ b.FrontId ].second; + exists &= this->m_CollisionSites[ b.FrontId ][ a.FrontId ].second; + if( !exists ) + { + this->m_CollisionSites[ a.FrontId ][ b.FrontId ].first = a.Vertex; + this->m_CollisionSites[ a.FrontId ][ b.FrontId ].second = true; + this->m_CollisionSites[ b.FrontId ][ a.FrontId ].first = b.Vertex; + this->m_CollisionSites[ b.FrontId ][ a.FrontId ].second = true; + + // Stop if one front is desired + if( this->m_StopAtOneFront ) + { + // Perform a depth-first iteration on front graph + unsigned long N = this->GetNumberOfSeeds( ); + unsigned long C = 0; + std::vector< bool > m( N, false ); + std::queue< unsigned long > q; + q.push( 0 ); + while( !q.empty( ) ) + { + unsigned long f = q.front( ); + q.pop( ); + + if( m[ f ] ) + continue; + m[ f ] = true; + C++; + + for( unsigned int n = 0; n < N; ++n ) + if( this->m_CollisionSites[ f ][ n ].second && !m[ n ] ) + q.push( n ); + + } // elihw + ret = ( C == N ); + + } // fi + + } // fi + + } // fi + return( ret ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +bool fpa::Base::Algorithm< T, B >:: +_CheckStopCondition( ) +{ + return( false ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +bool fpa::Base::Algorithm< T, B >:: +_UpdateResult( _TNode& n ) +{ + return( true ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +_InitializeMarks( ) +{ + this->m_Marks.clear( ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +bool fpa::Base::Algorithm< T, B >:: +_IsMarked( const TVertex& v ) const +{ + return( this->m_Marks.find( v ) != this->m_Marks.end( ) ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +typename fpa::Base::Algorithm< T, B >:: +_TFrontId fpa::Base::Algorithm< T, B >:: +_FrontId( const TVertex& v ) const +{ + typename _TMarks::const_iterator mIt = this->m_Marks.find( v ); + if( mIt != this->m_Marks.end( ) ) + return( mIt->second.FrontId ); + else + return( std::numeric_limits< _TFrontId >::max( ) ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +typename fpa::Base::Algorithm< T, B >:: +TVertex fpa::Base::Algorithm< T, B >:: +_Parent( const TVertex& v ) const +{ + typename _TMarks::const_iterator mIt = this->m_Marks.find( v ); + if( mIt == this->m_Marks.end( ) ) + { + TVertex v; + return( v ); + } + else + return( mIt->second.Parent ); +} + +// ------------------------------------------------------------------------- +template< class T, class B > +void fpa::Base::Algorithm< T, B >:: +_Mark( const _TNode& n ) +{ + this->m_Marks[ n.Vertex ] = n; +} + +#endif // __FPA__BASE__ALGORITHM__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Dijkstra.h b/lib/fpa/Base/Dijkstra.h new file mode 100644 index 0000000..a6212cf --- /dev/null +++ b/lib/fpa/Base/Dijkstra.h @@ -0,0 +1,120 @@ +#ifndef __FPA__BASE__DIJKSTRA__H__ +#define __FPA__BASE__DIJKSTRA__H__ + +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class V, class C, class VV, class VC > + class DijkstraTraits + { + public: + typedef V TVertex; + typedef C TResult; + typedef C TCost; + typedef VV TVertexValue; + typedef VC TVertexCmp; + typedef long TFrontId; + + class TNode + { + public: + TNode( ) + : Cost( 0 ) + { } + TNode( const TVertex& v, const TFrontId& f ) + : Vertex( v ), + Parent( v ), + Result( TResult( 0 ) ), + FrontId( f ), + Cost( TCost( 0 ) ) + { } + TNode( const TVertex& v, const TResult& r, const TFrontId& f ) + : Vertex( v ), + Parent( v ), + Result( r ), + FrontId( f ), + Cost( TCost( 0 ) ) + { } + virtual ~TNode( ) + { } + + // NOTE: stl::heaps work as maximum priority queues + bool operator<( const TNode& other ) const + { return( other.Cost < this->Cost ); } + + TVertex Vertex; + TVertex Parent; + TResult Result; + TFrontId FrontId; + TCost Cost; + }; + + typedef std::vector< TNode > TNodes; + }; + + /** + * Dijkstra is a front propagation algorithm that minimizes costs + */ + template< class V, class C, class VV, class VC, class B > + class Dijkstra + : public Algorithm< DijkstraTraits< V, C, VV, VC >, B > + { + public: + // Templated types + typedef V TVertex; + typedef C TCost; + typedef VV TVertexValue; + typedef B TBaseFilter; + typedef DijkstraTraits< V, C, VV, VC > TTraits; + + // Standard class typdedefs + typedef Dijkstra Self; + typedef Algorithm< TTraits, B > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + protected: + typedef typename TTraits::TFrontId _TFrontId; + typedef typename TTraits::TNode _TNode; + typedef typename TTraits::TNodes _TNodes; + + typedef std::vector< _TNode > _TQueue; + + public: + itkTypeMacro( Dijkstra, Base ); + + protected: + Dijkstra( ); + virtual ~Dijkstra( ); + + virtual void _InitializeQueue ( ); + virtual bool _IsQueueEmpty ( ) const; + virtual void _QueuePush ( const _TNode& n ); + virtual _TNode _QueuePop ( ); + virtual void _QueueClear ( ); + virtual bool _UpdateNeigh ( _TNode& nn, const _TNode& n ); + + private: + // Purposely not implemented + Dijkstra( const Self& ); + void operator=( const Self& ); + + private: + _TQueue m_Queue; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__BASE__DIJKSTRA__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Dijkstra.hxx b/lib/fpa/Base/Dijkstra.hxx new file mode 100644 index 0000000..4920ab5 --- /dev/null +++ b/lib/fpa/Base/Dijkstra.hxx @@ -0,0 +1,95 @@ +#ifndef __FPA__BASE__DIJKSTRA__HXX__ +#define __FPA__BASE__DIJKSTRA__HXX__ + +#include + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +fpa::Base::Dijkstra< V, C, VV, VC, B >:: +Dijkstra( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +fpa::Base::Dijkstra< V, C, VV, VC, B >:: +~Dijkstra( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +void fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_InitializeQueue( ) +{ + for( + typename _TNodes::const_iterator vIt = this->m_Seeds.begin( ); + vIt != this->m_Seeds.end( ); + vIt++ + ) + this->_QueuePush( *vIt ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +bool fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_IsQueueEmpty( ) const +{ + return( this->m_Queue.empty( ) ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +void fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_QueuePush( const _TNode& n ) +{ + this->m_Queue.push_back( n ); + std::push_heap( this->m_Queue.begin( ), this->m_Queue.end( ) ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +typename fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_TNode fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_QueuePop( ) +{ + _TNode n; + if( !( this->m_Queue.empty( ) ) ) + { + // n = *( this->m_Queue.begin( ) ); + n = this->m_Queue.front( ); + std::pop_heap( this->m_Queue.begin( ), this->m_Queue.end( ) ); + this->m_Queue.pop_back( ); + + } // fi + return( n ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +void fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_QueueClear( ) +{ + this->m_Queue.clear( ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +bool fpa::Base::Dijkstra< V, C, VV, VC, B >:: +_UpdateNeigh( _TNode& nn, const _TNode& n ) +{ + TCost nc = this->_Cost( nn.Vertex, n.Vertex ); + if( C( 0 ) <= nc ) + { + nn.Cost = n.Cost + nc; + nn.Result = nn.Cost; + return( true ); + } + else + return( false ); +} + +#endif // __FPA__BASE__DIJKSTRA__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Events.h b/lib/fpa/Base/Events.h new file mode 100644 index 0000000..b7b9f4e --- /dev/null +++ b/lib/fpa/Base/Events.h @@ -0,0 +1,136 @@ +#ifndef __FPA__BASE__EVENTS__H__ +#define __FPA__BASE__EVENTS__H__ + +#include +#include + +namespace fpa +{ + namespace Base + { + /** + * Evolution event. An event is generated when a vertex changes its + * state. + */ + template< class N > + class BaseEvent + : public itk::AnyEvent + { + public: + BaseEvent( ) + : itk::AnyEvent( ) + { } + BaseEvent( const N& n ) + : itk::AnyEvent( ), + Node( n ) + { } + virtual ~BaseEvent( ) + { } + + const char* GetEventName( ) const + { return( "fpa::Base::BaseEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const BaseEvent< N >* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new BaseEvent< N >( ) ); } + + public: + N Node; + }; + + /** + */ + template< class N > + class FrontEvent + : public BaseEvent< N > + { + public: + FrontEvent( ) + : BaseEvent< N >( ) + { } + FrontEvent( const N& n ) + : BaseEvent< N >( n ) + { } + virtual ~FrontEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::FrontEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const FrontEvent< N >* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new FrontEvent< N >( ) ); } + }; + + /** + */ + template< class N > + class MarkEvent + : public BaseEvent< N > + { + public: + MarkEvent( ) + : BaseEvent< N >( ) + { } + MarkEvent( const N& n ) + : BaseEvent< N >( n ) + { } + virtual ~MarkEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::MarkEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const MarkEvent< N >* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new MarkEvent< N >( ) ); } + }; + + /** + */ + template< class N > + class CollisionEvent + : public BaseEvent< N > + { + public: + CollisionEvent( ) + : BaseEvent< N >( ) + { } + CollisionEvent( const N& n ) + : BaseEvent< N >( n ) + { } + virtual ~CollisionEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::CollisionEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const CollisionEvent< N >* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new CollisionEvent< N >( ) ); } + }; + + /** + */ + template< class N > + class EndEvent + : public BaseEvent< N > + { + public: + EndEvent( ) + : BaseEvent< N >( ) + { } + virtual ~EndEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::EndEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const EndEvent< N >* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new EndEvent< N >( ) ); } + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__BASE__EVENTS__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/FastMarching.h b/lib/fpa/Base/FastMarching.h new file mode 100644 index 0000000..119d0f3 --- /dev/null +++ b/lib/fpa/Base/FastMarching.h @@ -0,0 +1,76 @@ +#ifndef __FPA__BASE__FASTMARCHING__H__ +#define __FPA__BASE__FASTMARCHING__H__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class V, class C, class VV, class VC, class B > + class FastMarching + : public Dijkstra< V, C, VV, VC, B > + { + public: + // Templated types + typedef V TVertex; + typedef C TCost; + typedef VV TVertexValue; + typedef B TBaseFilter; + + // Standard class typdedefs + typedef FastMarching Self; + typedef Dijkstra< V, C, VV, VC, B > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TTraits TTraits; + + protected: + typedef typename Superclass::_TFrontId _TFrontId; + typedef typename Superclass::_TNode _TNode; + typedef typename Superclass::_TNodes _TNodes; + struct _TNodeEikonalSort + { + bool operator()( const _TNode& a, const _TNode& b ) const + { return( a.Result < b.Result ); } + }; + + public: + itkTypeMacro( FastMarching, Dijkstra ); + + protected: + FastMarching( ); + virtual ~FastMarching( ); + + virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n ); + + /** + * Eikonal solution method + * + * This method returns the computed time and a flag indicating if the + * determinant was valid. + */ + bool _SolveEikonal( double& re, const _TNode& n, const double& F ); + + private: + // Purposely not implemented + FastMarching( const Self& ); + void operator=( const Self& ); + + protected: + static const double LargestValue; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__BASE__FASTMARCHING__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/FastMarching.hxx b/lib/fpa/Base/FastMarching.hxx new file mode 100644 index 0000000..ef1f634 --- /dev/null +++ b/lib/fpa/Base/FastMarching.hxx @@ -0,0 +1,121 @@ +#ifndef __FPA__BASE__FASTMARCHING__HXX__ +#define __FPA__BASE__FASTMARCHING__HXX__ + +#include +#include + +#include + +// --------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +const double fpa::Base::FastMarching< V, C, VV, VC, B >:: +LargestValue = std::numeric_limits< double >::max( ); + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +fpa::Base::FastMarching< V, C, VV, VC, B >:: +FastMarching( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +fpa::Base::FastMarching< V, C, VV, VC, B >:: +~FastMarching( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +bool fpa::Base::FastMarching< V, C, VV, VC, B >:: +_UpdateNeigh( _TNode& nn, const _TNode& n ) +{ + bool r = true; + if( !( this->_IsMarked( nn.Vertex ) ) ) + { + double re; + double F = double( this->_Value( nn.Vertex ) ); + r = this->_SolveEikonal( re, nn, F ); + if( r ) + { + nn.Result = TResult( re ); + nn.Cost = TCost( re ); + + } // fi + + } // fi + return( r ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class VV, class VC, class B > +bool fpa::Base::FastMarching< V, C, VV, VC, B >:: +_SolveEikonal( double& re, const _TNode& n, const double& F ) +{ + _TNodes neighs; + this->_Neighs( n, neighs ); + + // Compute interest nodes (i.e. the smallest valued in each dimension) + _TNodes minNodes; + unsigned int i = 0; + while( i < neighs.size( ) ) + { + V nm = neighs[ i++ ].Vertex; + V np = neighs[ i++ ].Vertex; + bool mm = this->_IsMarked( nm ); + bool mp = this->_IsMarked( np ); + TResult rm = this->_Result( nm ); + TResult rp = this->_Result( np ); + + _TNode Ti; + if( mm && mp ) + minNodes.push_back( + ( rm < rp )? + _TNode( nm, rm, n.FrontId ): + _TNode( np, rp, n.FrontId ) + ); + else if( !mm && mp ) + minNodes.push_back( _TNode( np, rp, n.FrontId ) ); + else if( mm && !mp ) + minNodes.push_back( _TNode( nm, rm, n.FrontId ) ); + + } // elihw + std::sort( minNodes.begin( ), minNodes.end( ), _TNodeEikonalSort( ) ); + + // Compute solution + double al = double( 0 ); + double be = double( 0 ); + double ga = double( -1 ) / ( F * F ); + re = Self::LargestValue; + + typename _TNodes::const_iterator nIt = minNodes.begin( ); + bool valid = true; + while( nIt != minNodes.end( ) && valid ) + { + if( !( re < double( nIt->Result ) ) ) + { + double d = this->_Norm( n.Vertex, nIt->Vertex ); + double s = double( 1 ) / ( d * d ); + double v = double( nIt->Result ); + + al += s; + be += v * s; + ga += v * v * s; + + double de = ( be * be ) - ( al * ga ); + valid = !( de < double( 0 ) ); + if( valid ) + re = ( be + std::sqrt( de ) ) / al; + nIt++; + } + else + nIt = minNodes.end( ); + + } // elihw + return( valid ); +} + +#endif // __FPA__BASE__FASTMARCHING__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Functors/InvertCostFunction.h b/lib/fpa/Base/Functors/InvertCostFunction.h new file mode 100644 index 0000000..faf4fb2 --- /dev/null +++ b/lib/fpa/Base/Functors/InvertCostFunction.h @@ -0,0 +1,56 @@ +#ifndef __FPA__BASE__FUNCTORS__INVERTCOSTFUNCTION__H__ +#define __FPA__BASE__FUNCTORS__INVERTCOSTFUNCTION__H__ + +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class C > + class InvertCostFunction + : public itk::FunctionBase< C, C > + { + public: + /// Type-related and pointers + typedef InvertCostFunction Self; + typedef itk::FunctionBase< C, C > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( InvertCostFunction, itkFunctionBase ); + + public: + virtual C Evaluate( const C& input ) const + { + return( C( 1 ) / ( C( 1 ) + C( input ) ) ); + } + + protected: + InvertCostFunction( ) + : Superclass( ) + { } + virtual ~InvertCostFunction( ) + { } + + private: + // Purposely not implemented + InvertCostFunction( const Self& ); + void operator=( const Self& ); + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __FPA__BASE__FUNCTORS__INVERTCOSTFUNCTION__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/RegionGrow.h b/lib/fpa/Base/RegionGrow.h new file mode 100644 index 0000000..e52fb62 --- /dev/null +++ b/lib/fpa/Base/RegionGrow.h @@ -0,0 +1,117 @@ +#ifndef __FPA__BASE__REGIONGROW__H__ +#define __FPA__BASE__REGIONGROW__H__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class V, class R, class VV, class VC > + class RegionGrowTraits + { + public: + typedef R TResult; + typedef V TVertex; + typedef VV TVertexValue; + typedef VC TVertexCmp; + + typedef bool TCost; + typedef long TFrontId; + + class TNode + { + public: + TNode( ) + { } + TNode( const TVertex& v, const TFrontId& f ) + : Vertex( v ), + Parent( v ), + FrontId( f ) + { } + TNode( const TVertex& v, const TResult& r, const TFrontId& f ) + : Vertex( v ), + Parent( v ), + Result( r ), + FrontId( f ) + { } + virtual ~TNode( ) + { } + + TVertex Vertex; + TVertex Parent; + TResult Result; + TFrontId FrontId; + }; + + typedef std::vector< TNode > TNodes; + }; + + /** + * Region grow is a front propagation with no costs. + */ + template< class V, class R, class VV, class VC, class B > + class RegionGrow + : public Algorithm< RegionGrowTraits< V, R, VV, VC >, B > + { + public: + typedef V TVertex; + typedef R TResult; + typedef VV TVertexValue; + typedef B TBaseFilter; + + /// Standard class typdedefs + typedef RegionGrowTraits< V, R, VV, VC > TTraits; + typedef RegionGrow Self; + typedef Algorithm< TTraits, B > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename TTraits::TCost TCost; + + protected: + typedef typename TTraits::TFrontId _TFrontId; + typedef typename TTraits::TNode _TNode; + typedef typename TTraits::TNodes _TNodes; + + private: + typedef std::queue< _TNode > _TQueue; + + public: + itkTypeMacro( RegionGrow, Base ); + + protected: + RegionGrow( ); + virtual ~RegionGrow( ); + + virtual bool _UpdateResult ( _TNode& n ); + virtual void _InitializeQueue ( ); + virtual bool _IsQueueEmpty ( ) const; + virtual void _QueuePush ( const _TNode& n ); + virtual _TNode _QueuePop ( ); + virtual void _QueueClear ( ); + virtual bool _UpdateNeigh ( _TNode& nn, const _TNode& n ); + + virtual bool _CheckMembership( const _TNode& n ) const = 0; + + private: + RegionGrow( const Self& ); // Not impl. + void operator=( const Self& ); // Not impl. + + private: + _TQueue m_Queue; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__BASE__REGIONGROW__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/RegionGrow.hxx b/lib/fpa/Base/RegionGrow.hxx new file mode 100644 index 0000000..bdc1ed3 --- /dev/null +++ b/lib/fpa/Base/RegionGrow.hxx @@ -0,0 +1,87 @@ +#ifndef __FPA__BASE__REGIONGROWING__HXX__ +#define __FPA__BASE__REGIONGROWING__HXX__ + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +fpa::Base::RegionGrow< V, R, VV, VC, B >:: +RegionGrow( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +fpa::Base::RegionGrow< V, R, VV, VC, B >:: +~RegionGrow( ) +{ +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +bool fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_UpdateResult( _TNode& n ) +{ + n.Result = R( this->_CheckMembership( n ) ); + return( n.Result ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +void fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_InitializeQueue( ) +{ + for( + typename _TNodes::const_iterator vIt = this->m_Seeds.begin( ); + vIt != this->m_Seeds.end( ); + vIt++ + ) + this->_QueuePush( *vIt ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +bool fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_IsQueueEmpty( ) const +{ + return( this->m_Queue.empty( ) ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +void fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_QueuePush( const _TNode& n ) +{ + this->m_Queue.push( n ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +typename fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_TNode fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_QueuePop( ) +{ + _TNode n = this->m_Queue.front( ); + this->m_Queue.pop( ); + return( n ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +void fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_QueueClear( ) +{ + while( this->m_Queue.size( ) > 0 ) + this->m_Queue.pop( ); +} + +// ------------------------------------------------------------------------- +template< class V, class R, class VV, class VC, class B > +bool fpa::Base::RegionGrow< V, R, VV, VC, B >:: +_UpdateNeigh( _TNode& nn, const _TNode& n ) +{ + return( true ); +} + +#endif // __FPA__BASE__REGIONGROWING__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/TreeExtractor.h b/lib/fpa/Base/TreeExtractor.h new file mode 100644 index 0000000..a0222a6 --- /dev/null +++ b/lib/fpa/Base/TreeExtractor.h @@ -0,0 +1,84 @@ +#ifndef __FPA__BASE__TREEEXTRACTOR__H__ +#define __FPA__BASE__TREEEXTRACTOR__H__ + +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class A > + class TreeExtractor + : public A + { + public: + /// Standard class typdedefs + typedef TreeExtractor Self; + typedef A Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef A TAlgorithm; + + typedef typename A::TCost TCost; + typedef typename A::TVertex TVertex; + typedef std::vector< TVertex > TVertices; + typedef typename A::TTraits::TVertexCmp TVertexCmp; + typedef std::map< TVertex, TVertex, TVertexCmp > TTree; + + protected: + typedef typename A::_TFrontId _TFrontId; + typedef typename A::_TNode _TNode; + typedef typename A::_TCollisionSitesRow _TCollisionSitesRow; + typedef typename A::_TCollisionSites _TCollisionSites; + typedef std::vector< unsigned long > _TRow; + typedef std::vector< _TRow > _TMatrix; + + public: + itkNewMacro( Self ); + itkTypeMacro( TreeExtractor, TAlgorithm ); + + itkGetConstMacro( Root, TVertex ); + itkGetConstMacro( Tree, TTree ); + itkSetMacro( Root, TVertex ); + + public: + unsigned int AddLeaf( const TVertex& leaf ); + unsigned int GetNumberOfLeafs( ) const; + const TVertex& GetLeaf( const unsigned int& id ) const; + + protected: + TreeExtractor( ); + virtual ~TreeExtractor( ); + + virtual void _BeforeLoop( ); + virtual void _AfterLoop( ); + + TVertices _Path( const TVertex& s, const TVertex& e ); + + private: + // Purporsely not implemented + TreeExtractor( const Self& ); + void operator=( const Self& ); + + protected: + TVertex m_Root; + TVertices m_Leafs; + TTree m_Tree; + + _TMatrix m_FrontPaths; + static const unsigned long INF_VALUE; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__BASE__TREEEXTRACTOR__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/TreeExtractor.hxx b/lib/fpa/Base/TreeExtractor.hxx new file mode 100644 index 0000000..cfeb17f --- /dev/null +++ b/lib/fpa/Base/TreeExtractor.hxx @@ -0,0 +1,250 @@ +#ifndef __FPA__BASE__TREEEXTRACTOR__HXX__ +#define __FPA__BASE__TREEEXTRACTOR__HXX__ + +#include + +// ------------------------------------------------------------------------- +template< class A > +const unsigned long fpa::Base::TreeExtractor< A >::INF_VALUE = + std::numeric_limits< unsigned long >::max( ) >> 1; + +// ------------------------------------------------------------------------- +template< class A > +unsigned int fpa::Base::TreeExtractor< A >:: +AddLeaf( const TVertex& leaf ) +{ + this->m_Leafs.push_back( leaf ); + this->Modified( ); + return( this->m_Leafs.size( ) - 1 ); +} + +// ------------------------------------------------------------------------- +template< class A > +unsigned int fpa::Base::TreeExtractor< A >:: +GetNumberOfLeafs( ) const +{ + return( this->m_Leafs.size( ) ); +} + +// ------------------------------------------------------------------------- +template< class A > +const typename fpa::Base::TreeExtractor< A >:: +TVertex& fpa::Base::TreeExtractor< A >:: +GetLeaf( const unsigned int& id ) const +{ + static const TVertex zero; + if( id < this->m_Leafs.size( ) ) + return( this->m_Leafs[ id ] ); + else + return( zero ); +} + +// ------------------------------------------------------------------------- +template< class A > +fpa::Base::TreeExtractor< A >:: +TreeExtractor( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class A > +fpa::Base::TreeExtractor< A >:: +~TreeExtractor( ) +{ +} + +// ------------------------------------------------------------------------- +template< class A > +void fpa::Base::TreeExtractor< A >:: +_BeforeLoop( ) +{ + this->Superclass::_BeforeLoop( ); + this->m_Tree.clear( ); +} + +// ------------------------------------------------------------------------- +template< class A > +void fpa::Base::TreeExtractor< A >:: +_AfterLoop( ) +{ + this->Superclass::_AfterLoop( ); + + // Prepare fronts graph using Floyd-Warshall + unsigned long nSeeds = this->GetNumberOfSeeds( ); + _TMatrix dist( nSeeds, _TRow( nSeeds, Self::INF_VALUE ) ); + this->m_FrontPaths = dist; + + for( unsigned long i = 0; i < nSeeds; ++i ) + { + for( unsigned long j = 0; j < nSeeds; ++j ) + { + if( this->m_CollisionSites[ i ][ j ].second ) + { + dist[ i ][ j ] = 1; + dist[ j ][ i ] = 1; + this->m_FrontPaths[ i ][ j ] = j; + this->m_FrontPaths[ j ][ i ] = i; + + } // fi + + } // rof + dist[ i ][ i ] = 0; + this->m_FrontPaths[ i ][ i ] = i; + + } // rof + for( unsigned long k = 0; k < nSeeds; ++k ) + { + for( unsigned long i = 0; i < nSeeds; ++i ) + { + for( unsigned long j = 0; j < nSeeds; ++j ) + { + if( ( dist[ i ][ k ] + dist[ k ][ j ] ) < dist[ i ][ j ] ) + { + dist[ i ][ j ] = dist[ i ][ k ] + dist[ k ][ j ]; + this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ]; + + } // fi + + } // rof + + } // rof + + } // rof + + // Build tree + this->m_Tree.clear( ); + typename TVertices::const_iterator vIt = this->m_Leafs.begin( ); + for( ; vIt != this->m_Leafs.end( ); ++vIt ) + { + TVertices path = this->_Path( this->m_Root, *vIt ); + typename TVertices::const_iterator pIt = path.begin( ); + TVertex parent = *pIt; + for( ; pIt != path.end( ); ++pIt ) + { + this->m_Tree[ *pIt ] = parent; + parent = *pIt; + + } // rof + + } // rof +} + +// ------------------------------------------------------------------------- +template< class A > +typename fpa::Base::TreeExtractor< A >:: +TVertices fpa::Base::TreeExtractor< A >:: +_Path( const TVertex& s, const TVertex& e ) +{ + TVertices pS, pE; + + // Check if both vertices belong to the same front or if collisions + // should be analyzed. + _TFrontId fs = this->_FrontId( s ); + _TFrontId fe = this->_FrontId( e ); + if( fs == fe ) + { + // Path from start vertex until seed + TVertex vIt = s; + TVertex parent = this->_Parent( vIt ); + while( parent != vIt ) + { + pS.push_back( vIt ); + vIt = parent; + parent = this->_Parent( vIt ); + + } // elihw + pS.push_back( vIt ); + + // Path from end vertex until seed + vIt = e; + parent = this->_Parent( vIt ); + while( parent != vIt ) + { + pE.push_back( vIt ); + vIt = parent; + parent = this->_Parent( vIt ); + + } // elihw + pE.push_back( vIt ); + + // Erase duplicate vertices + if( pS.size( ) > 0 && pE.size( ) > 0 ) + { + bool cont = ( pS.back( ) == pE.back( ) ); + bool enter = false; + TVertex lastVertex = pS.back( ); + while( cont ) + { + enter = true; + + lastVertex = pS.back( ); + pS.pop_back( ); + pE.pop_back( ); + + cont = ( pS.size( ) > 0 && pE.size( ) > 0 ); + if( cont ) + cont = ( pS.back( ) == pE.back( ) ); + + } // elihw + if( enter ) + pS.push_back( lastVertex ); + + } // fi + + // Contatenate both paths + pS.insert( pS.end( ), pE.rbegin( ), pE.rend( ) ); + } + else + { + // Use this->m_FrontPaths from Floyd-Warshall + if( this->m_FrontPaths[ fs ][ fe ] != Self::INF_VALUE ) + { + // Compute front path + std::vector< unsigned long > fpath; + fpath.push_back( fs ); + while( fs != fe ) + { + fs = this->m_FrontPaths[ fs ][ fe ]; + fpath.push_back( fs ); + + } // elihw + + // Continue only if both fronts are connected + unsigned int N = fpath.size( ); + if( 0 < N ) + { + // First path: from start vertex to first collision + pS = this->_Path( + s, this->m_CollisionSites[ fpath[ 0 ] ][ fpath[ 1 ] ].first + ); + + // Intermediary paths + for( unsigned int i = 1; i < N - 1; ++i ) + { + pE = this->_Path( + this->m_CollisionSites[ fpath[ i ] ][ fpath[ i - 1 ] ].first, + this->m_CollisionSites[ fpath[ i ] ][ fpath[ i + 1 ] ].first + ); + pS.insert( pS.end( ), pE.begin( ), pE.end( ) ); + + } // rof + + // Final path: from last collision to end point + pE = this->_Path( + this->m_CollisionSites[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, + e + ); + pS.insert( pS.end( ), pE.begin( ), pE.end( ) ); + + } // fi + + } // fi + + } // fi + return( pS ); +} + +#endif // __FPA__BASE__TREEEXTRACTOR__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Common.cxx.in b/lib/fpa/Common.cxx.in new file mode 100644 index 0000000..2b1a193 --- /dev/null +++ b/lib/fpa/Common.cxx.in @@ -0,0 +1,9 @@ +#include + +// ------------------------------------------------------------------------- +std::string fpa::GetVersion( ) +{ + return( "@VERSION_STRING@" ); +} + +// eof - $RCSfile$ diff --git a/lib/fpa/Common.h b/lib/fpa/Common.h new file mode 100644 index 0000000..c5bf168 --- /dev/null +++ b/lib/fpa/Common.h @@ -0,0 +1,15 @@ +#ifndef __FPA__COMMON__H__ +#define __FPA__COMMON__H__ + +#include +#include + +namespace fpa +{ + std::string FrontAlgorithms_EXPORT GetVersion( ); + +} // ecapseman + +#endif // __FPA__COMMON__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Algorithm.h b/lib/fpa/Image/Algorithm.h new file mode 100644 index 0000000..289eb6c --- /dev/null +++ b/lib/fpa/Image/Algorithm.h @@ -0,0 +1,105 @@ +#ifndef __FPA__IMAGE__ALGORITHM__H__ +#define __FPA__IMAGE__ALGORITHM__H__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * A generic front propagation algorithm were vertices are image pixels. + * + * @param I Input image type + * @param A Base algorithm (RegionGrow, Dijkstra or FastMarching) + */ + template< class I, class A > + class Algorithm + : public A + { + public: + /// Standard class typdedefs + typedef Algorithm Self; + typedef A Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /// Template input values + typedef I TInputImage; + typedef A TBaseAlgorithm; + + typedef typename A::TTraits TTraits; + typedef typename TTraits::TCost TCost; + typedef typename TTraits::TResult TResult; + typedef typename TTraits::TVertex TVertex; + typedef typename TTraits::TVertexValue TVertexValue; + + typedef itk::Image< TResult, I::ImageDimension > TOutputImage; + typedef itk::FunctionBase< TCost, TCost > TCostConversionFunction; + + protected: + typedef typename TTraits::TFrontId _TFrontId; + typedef typename TTraits::TNode _TNode; + typedef typename TTraits::TNodes _TNodes; + + private: + typedef itk::Image< bool, I::ImageDimension > _TMarks; + typedef itk::Image< _TFrontId, I::ImageDimension > _TFrontsIds; + typedef itk::Image< TVertex, I::ImageDimension > _TParents; + + public: + itkTypeMacro( Algorithm, TAlgorithm ); + + /// Set/Get + itkGetConstMacro( NeighborhoodOrder, unsigned int ); + itkGetConstObjectMacro( CostConversion, TCostConversionFunction ); + itkGetObjectMacro( CostConversion, TCostConversionFunction ); + + itkSetMacro( NeighborhoodOrder, unsigned int ); + itkSetObjectMacro( CostConversion, TCostConversionFunction ); + + protected: + Algorithm( ); + virtual ~Algorithm( ); + + /// Base interface + virtual bool _UpdateResult( _TNode& n ); + + /// Pure virtual interface: vertices + virtual unsigned long _NumberOfVertices ( ) const; + virtual TVertexValue _Value ( const TVertex& v ) const; + virtual TResult _Result ( const TVertex& v ) const; + + /// Pure virtual interface: edges + virtual double _Norm ( const TVertex& a, const TVertex& b ) const; + virtual bool _Edge ( const TVertex& a, const TVertex& b ) const; + virtual TCost _Cost ( const TVertex& a, const TVertex& b ) const; + + /// Pure virtual interface: neighborhood + virtual void _Neighs ( const _TNode& n, _TNodes& N ) const; + virtual void _NeighsInDim ( const _TNode& n, + const unsigned int& d, + _TNodes& N ); + + /// Pure virtual interface: results + virtual void _InitializeResults ( ); + + private: + Algorithm( const Self& ); // Not impl. + void operator=( const Self& ); // Not impl. + + protected: + unsigned int m_NeighborhoodOrder; + typename TCostConversionFunction::Pointer m_CostConversion; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IMAGE__ALGORITHM__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Algorithm.hxx b/lib/fpa/Image/Algorithm.hxx new file mode 100644 index 0000000..7800a54 --- /dev/null +++ b/lib/fpa/Image/Algorithm.hxx @@ -0,0 +1,180 @@ +#ifndef __FPA__IMAGE__ALGORITHM__HXX__ +#define __FPA__IMAGE__ALGORITHM__HXX__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class I, class A > +fpa::Image::Algorithm< I, A >:: +Algorithm( ) + : Superclass( ), + m_NeighborhoodOrder( 1 ) +{ +} + +// ------------------------------------------------------------------------- +template< class I, class A > +fpa::Image::Algorithm< I, A >:: +~Algorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template< class I, class A > +bool fpa::Image::Algorithm< I, A >:: +_UpdateResult( _TNode& n ) +{ + bool ret = this->Superclass::_UpdateResult( n ); + this->GetOutput( )->SetPixel( n.Vertex, n.Result ); + return( ret ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +unsigned long fpa::Image::Algorithm< I, A >:: +_NumberOfVertices( ) const +{ + return( + this->GetInput( )->GetLargestPossibleRegion( ).GetNumberOfPixels( ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +typename fpa::Image::Algorithm< I, A >:: +TVertexValue fpa::Image::Algorithm< I, A >:: +_Value( const TVertex& v ) const +{ + return( this->GetInput( )->GetPixel( v ) ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +typename fpa::Image::Algorithm< I, A >:: +TResult fpa::Image::Algorithm< I, A >:: +_Result( const TVertex& v ) const +{ + return( this->GetOutput( )->GetPixel( v ) ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +double fpa::Image::Algorithm< I, A >:: +_Norm( const TVertex& a, const TVertex& b ) const +{ + typename I::PointType pa, pb; + this->GetInput( )->TransformIndexToPhysicalPoint( a, pa ); + this->GetInput( )->TransformIndexToPhysicalPoint( b, pb ); + return( double( pa.EuclideanDistanceTo( pb ) ) ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +bool fpa::Image::Algorithm< I, A >:: +_Edge( const TVertex& a, const TVertex& b ) const +{ + unsigned long dist = 0; + for( unsigned int d = 0; d < I::ImageDimension; d++ ) + dist += std::abs( long( a[ d ] ) - long( b[ d ] ) ); + if( this->m_NeighborhoodOrder == 1 ) + return( dist == 0 || dist == 1 ); + else + return( dist == 0 || dist == 1 || dist == 2 ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +typename fpa::Image::Algorithm< I, A >:: +TCost fpa::Image::Algorithm< I, A >:: +_Cost( const TVertex& a, const TVertex& b ) const +{ + static const TCost INF_COST = std::numeric_limits< TCost >::max( ); + if( this->_Edge( a, b ) ) + { + TCost c = TCost( this->GetInput( )->GetPixel( b ) ); + if( this->m_CostConversion.IsNotNull( ) ) + return( this->m_CostConversion->Evaluate( c ) ); + else + return( c ); + } + else + return( INF_COST ); +} + +// ------------------------------------------------------------------------- +template< class I, class A > +void fpa::Image::Algorithm< I, A >:: +_Neighs( const _TNode& n, _TNodes& N ) const +{ + typename I::RegionType reg = this->GetInput( )->GetRequestedRegion( ); + + N.clear( ); + if( this->m_NeighborhoodOrder == 1 ) + { + for( unsigned int d = 0; d < I::ImageDimension; d++ ) + { + for( int i = -1; i <= 1; i += 2 ) + { + TVertex v = n.Vertex; + v[ d ] += i; + if( reg.IsInside( v ) ) + N.push_back( _TNode( v, n.FrontId ) ); + else + N.push_back( n ); + + } // rof + + } // rof + } + else + { + typedef itk::ConstNeighborhoodIterator< I > TNeighIt; + typename I::SizeType nSize; + nSize.Fill( 1 ); + + TNeighIt nIt( nSize, this->GetInput( ), reg ); + nIt.SetLocation( n.Vertex ); + for( unsigned int i = 0; i < nIt.Size( ); i++ ) + { + TVertex idxN = nIt.GetIndex( i ); + if( idxN == n.Vertex ) + continue; + if( !reg.IsInside( idxN ) ) + continue; + N.push_back( _TNode( idxN, n.FrontId ) ); + + } // rof + + } // fi +} + +// ------------------------------------------------------------------------- +template< class I, class A > +void fpa::Image::Algorithm< I, A >:: +_NeighsInDim( const _TNode& n, const unsigned int& d, _TNodes& N ) +{ + typename I::RegionType reg = this->GetInput( )->GetRequestedRegion( ); + + N.clear( ); + for( int i = -1; i <= 1; i += 2 ) + { + TVertex v = n.Vertex; + v[ d ] += i; + if( reg.IsInside( v ) ) + N.push_back( _TNode( v, n.FrontId ) ); + + } // rof +} + +// ------------------------------------------------------------------------- +template< class I, class A > +void fpa::Image::Algorithm< I, A >:: +_InitializeResults( ) +{ +} + +#endif // __FPA__IMAGE__ALGORITHM__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Dijkstra.h b/lib/fpa/Image/Dijkstra.h new file mode 100644 index 0000000..1464d3c --- /dev/null +++ b/lib/fpa/Image/Dijkstra.h @@ -0,0 +1,58 @@ +#ifndef __FPA__IMAGE__DIJKSTRA__H__ +#define __FPA__IMAGE__DIJKSTRA__H__ + +#include +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * @param I Input image type + */ + template< class I, class C > + class Dijkstra + : public Algorithm< I, fpa::Base::Dijkstra< typename I::IndexType, C, typename I::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, itk::Image< C, I::ImageDimension > > > > + { + public: + // Standard class typdedefs + typedef typename I::IndexType TVertex; + typedef typename I::PixelType TVertexValue; + typedef itk::Image< C, I::ImageDimension > TCostImage; + typedef itk::ImageToImageFilter< I, TCostImage > TBaseFilter; + typedef fpa::Base::Dijkstra< TVertex, C, TVertexValue, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, TBaseFilter > TBaseAlgorithm; + + typedef Dijkstra Self; + typedef Algorithm< I, TBaseAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( Dijkstra, fpaBaseDijkstra ); + + protected: + Dijkstra( ) + : Superclass( ) + { } + virtual ~Dijkstra( ) + { } + + private: + // Purposely not implemented + Dijkstra( const Self& ); + void operator=( const Self& ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__IMAGE__DIJKSTRA__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/FastMarching.h b/lib/fpa/Image/FastMarching.h new file mode 100644 index 0000000..d86d99d --- /dev/null +++ b/lib/fpa/Image/FastMarching.h @@ -0,0 +1,58 @@ +#ifndef __FPA__IMAGE__FASTMARCHING__H__ +#define __FPA__IMAGE__FASTMARCHING__H__ + +#include +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * @param I Input image type + */ + template< class I, class C > + class FastMarching + : public Algorithm< I, fpa::Base::FastMarching< typename I::IndexType, C, typename I::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, itk::Image< C, I::ImageDimension > > > > + { + public: + // Standard class typdedefs + typedef typename I::IndexType TVertex; + typedef typename I::PixelType TVertexValue; + typedef itk::Image< C, I::ImageDimension > TCostImage; + typedef itk::ImageToImageFilter< I, TCostImage > TBaseFilter; + typedef fpa::Base::FastMarching< TVertex, C, TVertexValue, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, TBaseFilter > TBaseAlgorithm; + + typedef FastMarching Self; + typedef Algorithm< I, TBaseAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( FastMarching, fpaBaseFastMarching ); + + protected: + FastMarching( ) + : Superclass( ) + { } + virtual ~FastMarching( ) + { } + + private: + // Purposely not implemented + FastMarching( const Self& ); + void operator=( const Self& ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__IMAGE__FASTMARCHING__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Functors/RegionGrowAllBelongsFunction.h b/lib/fpa/Image/Functors/RegionGrowAllBelongsFunction.h new file mode 100644 index 0000000..ca52050 --- /dev/null +++ b/lib/fpa/Image/Functors/RegionGrowAllBelongsFunction.h @@ -0,0 +1,78 @@ +#ifndef __FPA__IMAGE__FUNCTORS__REGIONGROWALLBELONGSFUNCTION__H__ +#define __FPA__IMAGE__FUNCTORS__REGIONGROWALLBELONGSFUNCTION__H__ + +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class I > + class RegionGrowAllBelongsFunction + : public itk::ImageFunction< I, bool > + { + public: + /// Type-related and pointers + typedef RegionGrowAllBelongsFunction Self; + typedef itk::ImageFunction< I, bool > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + // Superclass' types + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::InputPixelType InputPixelType; + typedef typename Superclass::OutputType OutputType; + typedef typename Superclass::CoordRepType CoordRepType; + typedef typename Superclass::IndexType IndexType; + typedef typename Superclass::IndexValueType IndexValueType; + typedef typename Superclass::ContinuousIndexType ContinuousIndexType; + typedef typename Superclass::PointType PointType; + + public: + itkNewMacro( Self ); + itkTypeMacro( RegionGrowAllBelongsFunction, itkImageFunction ); + + public: + virtual OutputType Evaluate( const PointType& point ) const + { + IndexType index; + this->ConvertPointToNearestIndex( point, index ); + return( this->EvaluateAtIndex( index ) ); + } + virtual OutputType EvaluateAtContinuousIndex( + const ContinuousIndexType& cindex + ) const + { + IndexType index; + this->ConvertContinuousIndexToNearestIndex( cindex, index ); + return( this->EvaluateAtIndex( index ) ); + } + virtual OutputType EvaluateAtIndex( const IndexType& index ) const + { return( true ); } + + protected: + RegionGrowAllBelongsFunction( ) + : Superclass( ) + { } + virtual ~RegionGrowAllBelongsFunction( ) + { } + + private: + // Purposely not implemented + RegionGrowAllBelongsFunction( const Self& ); + void operator=( const Self& ); + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __FPA__IMAGE__FUNCTORS__REGIONGROWALLBELONGSFUNCTION__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Functors/RegionGrowThresholdFunction.h b/lib/fpa/Image/Functors/RegionGrowThresholdFunction.h new file mode 100644 index 0000000..00a813f --- /dev/null +++ b/lib/fpa/Image/Functors/RegionGrowThresholdFunction.h @@ -0,0 +1,96 @@ +#ifndef __FPA__IMAGE__FUNCTORS__REGIONGROWTHRESHOLDFUNCTION__H__ +#define __FPA__IMAGE__FUNCTORS__REGIONGROWTHRESHOLDFUNCTION__H__ + +#include +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class I > + class RegionGrowThresholdFunction + : public RegionGrowAllBelongsFunction< I > + { + public: + /// Type-related and pointers + typedef RegionGrowThresholdFunction Self; + typedef RegionGrowAllBelongsFunction< I > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + // Superclass' types + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::InputPixelType InputPixelType; + typedef typename Superclass::OutputType OutputType; + typedef typename Superclass::CoordRepType CoordRepType; + typedef typename Superclass::IndexType IndexType; + typedef typename Superclass::IndexValueType IndexValueType; + typedef typename Superclass::ContinuousIndexType ContinuousIndexType; + typedef typename Superclass::PointType PointType; + + public: + itkNewMacro( Self ); + itkTypeMacro( + RegionGrowThresholdFunction, + RegionGrowAllBelongsFunction + ); + + itkGetConstMacro( LowerThreshold, InputPixelType ); + itkGetConstMacro( UpperThreshold, InputPixelType ); + + itkSetMacro( LowerThreshold, InputPixelType ); + itkSetMacro( UpperThreshold, InputPixelType ); + + public: + virtual OutputType EvaluateAtIndex( const IndexType& index ) const + { + const I* img = this->GetInputImage( ); + if( img != NULL ) + { + if( this->IsInsideBuffer( index ) ) + { + InputPixelType v = img->GetPixel( index ); + return( + this->m_LowerThreshold <= v && + v <= this->m_UpperThreshold + ); + + } // fi + + } // fi + return( false ); + } + + protected: + RegionGrowThresholdFunction( ) + : Superclass( ), + m_LowerThreshold( std::numeric_limits< InputPixelType >::min( ) ), + m_UpperThreshold( std::numeric_limits< InputPixelType >::max( ) ) + { } + virtual ~RegionGrowThresholdFunction( ) + { } + + private: + // Purposely not implemented + RegionGrowThresholdFunction( const Self& ); + void operator=( const Self& ); + + protected: + InputPixelType m_LowerThreshold; + InputPixelType m_UpperThreshold; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __FPA__IMAGE__FUNCTORS__REGIONGROWTHRESHOLDFUNCTION__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/RegionGrow.h b/lib/fpa/Image/RegionGrow.h new file mode 100644 index 0000000..d65fa6c --- /dev/null +++ b/lib/fpa/Image/RegionGrow.h @@ -0,0 +1,83 @@ +#ifndef __FPA__IMAGE__REGIONGROW__H__ +#define __FPA__IMAGE__REGIONGROW__H__ + +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * @param I Input image type + */ + template< class I > + class RegionGrow + : public Algorithm< I, fpa::Base::RegionGrow< typename I::IndexType, typename I::PixelType, typename I::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, I > > > + { + public: + // Standard class typdedefs + typedef typename I::IndexType TVertex; + typedef typename I::PixelType TResult; + typedef typename I::PixelType TVertexValue; + typedef itk::ImageToImageFilter< I, I > TBaseFilter; + typedef fpa::Base::RegionGrow< TVertex, TResult, TVertexValue, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, TBaseFilter > TBaseAlgorithm; + + typedef RegionGrow Self; + typedef Algorithm< I, TBaseAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef itk::ImageFunction< I, bool > TMembershipFunction; + + public: + itkNewMacro( Self ); + itkTypeMacro( RegionGrow, fpaBaseRegionGrow ); + + itkGetObjectMacro( MembershipFunction, TMembershipFunction ); + itkSetObjectMacro( MembershipFunction, TMembershipFunction ); + + protected: + RegionGrow( ) + : Superclass( ), + m_MembershipFunction( NULL ) + { } + virtual ~RegionGrow( ) + { } + + virtual bool _CheckMembership( + const typename Superclass::_TNode& n + ) const + { + if( this->m_MembershipFunction.IsNotNull( ) ) + { + if( + this->m_MembershipFunction->GetInputImage( ) != + this->GetInput( ) + ) + this->m_MembershipFunction->SetInputImage( this->GetInput( ) ); + return( this->m_MembershipFunction->EvaluateAtIndex( n.Vertex ) ); + } + else + return( false ); + } + + private: + // Purposely not implemented + RegionGrow( const Self& ); + void operator=( const Self& ); + + protected: + typename TMembershipFunction::Pointer m_MembershipFunction; + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__IMAGE__REGIONGROW__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.h b/lib/fpa/Image/RegionGrowWithMultipleThresholds.h new file mode 100644 index 0000000..11066fe --- /dev/null +++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.h @@ -0,0 +1,70 @@ +#ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__H__ +#define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__H__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * @param I Input image type + */ + template< class I > + class RegionGrowWithMultipleThresholds + : public RegionGrow< I > + { + public: + typedef RegionGrowWithMultipleThresholds Self; + typedef RegionGrow< I > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename I::PixelType TPixel; + + typedef std::map< TPixel, unsigned long > THistogram; + typedef typename Superclass::TBaseAlgorithm TBaseAlgorithm; + + protected: + typedef typename TBaseAlgorithm::_TNode _TNode; + + public: + itkNewMacro( Self ); + itkTypeMacro( RegionGrowWithMultipleThresholds, RegionGrow ); + + itkGetConstMacro( DerivativeThreshold, double ); + itkSetMacro( DerivativeThreshold, double ); + + public: + void AddThreshold( const TPixel& v ); + void AddThresholds( + const TPixel& t0, const TPixel& t1, + const unsigned int& samples + ); + + protected: + RegionGrowWithMultipleThresholds( ); + virtual ~RegionGrowWithMultipleThresholds( ); + + virtual bool _UpdateResult( _TNode& n ); + virtual void _AfterLoop( ); + + private: + RegionGrowWithMultipleThresholds( const Self& ); // Not impl. + void operator=( const Self& ); // Not impl. + + protected: + double m_DerivativeThreshold; + THistogram m_Histogram; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx new file mode 100644 index 0000000..b4d4331 --- /dev/null +++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx @@ -0,0 +1,137 @@ +#ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__ +#define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__ + +#include + +// ------------------------------------------------------------------------- +template< class I > +void fpa::Image::RegionGrowWithMultipleThresholds< I >:: +AddThreshold( const TPixel& v ) +{ + this->m_Histogram[ v ] = 0; + + typedef + fpa::Image::Functors::RegionGrowThresholdFunction< I > + TFunction; + TFunction* function = + dynamic_cast< TFunction* >( this->GetMembershipFunction( ) ); + if( function != NULL ) + { + function->SetLowerThreshold( this->m_Histogram.begin( )->first ); + function->SetUpperThreshold( this->m_Histogram.rbegin( )->first ); + + } // fi + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class I > +void fpa::Image::RegionGrowWithMultipleThresholds< I >:: +AddThresholds( + const TPixel& t0, + const TPixel& t1, + const unsigned int& samples + ) +{ + double off = + ( double( t1 ) - double( t0 ) ) / ( double( samples ) - double( 1 ) ); + for( unsigned int s = 0; s < samples; s++ ) + this->AddThreshold( TPixel( ( double( s ) * off ) + double( t0 ) ) ); +} + +// ------------------------------------------------------------------------- +template< class I > +fpa::Image::RegionGrowWithMultipleThresholds< I >:: +RegionGrowWithMultipleThresholds( ) + : Superclass( ), + m_DerivativeThreshold( double( 3 ) ) +{ + typedef + fpa::Image::Functors::RegionGrowThresholdFunction< I > + TFunction; + typename TFunction::Pointer function = TFunction::New( ); + this->SetMembershipFunction( function ); +} + +// ------------------------------------------------------------------------- +template< class I > +fpa::Image::RegionGrowWithMultipleThresholds< I >:: +~RegionGrowWithMultipleThresholds( ) +{ +} + +// ------------------------------------------------------------------------- +template< class I > +bool fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_UpdateResult( _TNode& n ) +{ + bool ret = this->Superclass::_UpdateResult( n ); + + if( ret ) + { + TPixel v = TPixel( this->_Cost( n.Vertex, n.Vertex ) ); + + typename THistogram::reverse_iterator hIt = this->m_Histogram.rbegin( ); + while( hIt != this->m_Histogram.rend( ) ) + { + if( v <= hIt->first ) + { + hIt->second += 1; + hIt++; + } + else + hIt = this->m_Histogram.rend( ); + + } // elihw + this->GetOutput( )->SetPixel( n.Vertex, v ); + + } // fi + return( ret ); +} + +// ------------------------------------------------------------------------- +template< class I > +void fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_AfterLoop( ) +{ + typename THistogram::iterator prevIt = this->m_Histogram.begin( ); + typename THistogram::iterator currIt = prevIt; currIt++; + typename THistogram::iterator nextIt = currIt; nextIt++; + double prev_d1 = double( 0 ); + for( ; nextIt != this->m_Histogram.end( ); ++prevIt, ++currIt, ++nextIt ) + { + double d1 = double( nextIt->second ) - double( prevIt->second ); + d1 /= double( nextIt->first ) - double( prevIt->first ); + + std::cout + << currIt->first << " " + << currIt->second << " " + << d1 << " " + << ( d1 - prev_d1 ) << std::endl; + + prev_d1 = d1; + + } // rof + + + /* + double prev = double( 0 ); + while( hIt != this->m_Histogram.end( ) ) + { + double curr = double( hIt->second ); + double deri = curr - prev; + + if( hIt != this->m_Histogram.begin( ) ) + std::cout << hIt->first << " " << curr << " " << deri << std::endl; + + prev = curr; + hIt++; + + } // rof + */ + this->Superclass::_AfterLoop( ); +} + +#endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/Image2DObserver.h b/lib/fpa/VTK/Image2DObserver.h new file mode 100644 index 0000000..68dcae0 --- /dev/null +++ b/lib/fpa/VTK/Image2DObserver.h @@ -0,0 +1,91 @@ +#ifndef __FPA__VTK__IMAGE2DOBSERVER__H__ +#define __FPA__VTK__IMAGE2DOBSERVER__H__ + +#include +#include +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace VTK + { + /** + */ + template< class F, class R > + class Image2DObserver + : public itk::Command + { + public: + typedef Image2DObserver Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef F TFilter; + typedef R TRenderWindow; + typedef typename TFilter::TInputImage TImage; + typedef typename TFilter::TVertex TVertex; + + typedef std::set< TVertex > TVertices; + + public: + itkNewMacro( Self ); + itkTypeMacro( Image2DObserver, itkCommand ); + + public: + void SetImage( const TImage* img, R* rw ); + void SetPixel( + typename TImage::IndexType idx, + unsigned char red, + unsigned char green, + unsigned char blue, + unsigned char alpha + ); + void Render( ); + void Execute( itk::Object* c, const itk::EventObject& e ) + { this->Execute( ( const itk::Object* )( c ), e ); } + void Execute( const itk::Object* c, const itk::EventObject& e ); + + protected: + Image2DObserver( ) + : Superclass( ), + m_RenderWindow( NULL ), + m_Number( 0 ), + m_Percentage( 0 ) + { + this->m_Stencil = TImageData::New( ); + this->m_StencilActor = TImageActor::New( ); + } + virtual ~Image2DObserver( ) + { } + + private: + Image2DObserver( const Self& ); // Not impl. + void operator=( const Self& ); // Not impl. + + protected: + typedef vtkSmartPointer< vtkImageActor > TImageActor; + typedef vtkSmartPointer< vtkImageData > TImageData; + + typename TImage::ConstPointer m_Image; + + R* m_RenderWindow; + TImageData m_Stencil; + TImageActor m_StencilActor; + unsigned long m_Number; + unsigned long m_Percentage; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__VTK__IMAGE2DOBSERVER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/Image2DObserver.hxx b/lib/fpa/VTK/Image2DObserver.hxx new file mode 100644 index 0000000..7ebe526 --- /dev/null +++ b/lib/fpa/VTK/Image2DObserver.hxx @@ -0,0 +1,188 @@ +#ifndef __FPA__VTK__IMAGE2DOBSERVER__HXX__ +#define __FPA__VTK__IMAGE2DOBSERVER__HXX__ + +#include +#include +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image2DObserver< F, R >:: +SetImage( const TImage* img, R* rw ) +{ + this->m_Image = img; + this->m_RenderWindow = rw; + + unsigned int minD = TImage::ImageDimension; + minD = ( minD < 3 )? minD: 3; + + int e[ 6 ] = { 0 }; + typename TImage::RegionType reg = this->m_Image->GetRequestedRegion( ); + for( unsigned int i = 0; i < minD; i++ ) + { + e[ ( i << 1 ) + 0 ] = reg.GetIndex( )[ i ]; + e[ ( i << 1 ) + 1 ] = reg.GetIndex( )[ i ] + reg.GetSize( )[ i ] - 1; + + } // rof + + typename TImage::SpacingType spac = this->m_Image->GetSpacing( ); + double s[ 3 ] = { 1, 1, 1 }; + for( unsigned int i = 0; i < minD; i++ ) + s[ i ] = double( spac[ i ] ); + + typename TImage::PointType orig = this->m_Image->GetOrigin( ); + double o[ 3 ] = { 0 }; + for( unsigned int i = 0; i < minD; i++ ) + o[ i ] = double( orig[ i ] ); + + this->m_Stencil->SetExtent( e ); + this->m_Stencil->SetSpacing( s ); + this->m_Stencil->SetOrigin( o ); + this->m_Stencil->AllocateScalars( VTK_UNSIGNED_CHAR, 4 ); + for( unsigned int i = 0; i < 3; i++ ) + this->m_Stencil->GetPointData( )-> + GetScalars( )->FillComponent( i, 255 ); + this->m_Stencil->GetPointData( )->GetScalars( )->FillComponent( 3, 0 ); + + this->m_StencilActor->SetInputData( this->m_Stencil ); + this->m_StencilActor->InterpolateOff( ); + + double nPix = + double( this->m_Image->GetRequestedRegion( ).GetNumberOfPixels( ) ); + this->m_Percentage = ( unsigned long )( nPix * 0.01 ); + this->m_Number = 0; + + if( this->m_RenderWindow != NULL ) + { + vtkRenderer* ren = + this->m_RenderWindow->GetRenderers( )->GetFirstRenderer( ); + ren->AddActor( this->m_StencilActor ); + this->Render( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image2DObserver< F, R >:: +SetPixel( + typename TImage::IndexType idx, + unsigned char red, + unsigned char green, + unsigned char blue, + unsigned char alpha + ) +{ + this->m_Stencil->SetScalarComponentFromDouble + ( idx[ 0 ], idx[ 1 ], 0, 0, red ); + this->m_Stencil->SetScalarComponentFromDouble + ( idx[ 0 ], idx[ 1 ], 0, 1, green ); + this->m_Stencil->SetScalarComponentFromDouble + ( idx[ 0 ], idx[ 1 ], 0, 2, blue ); + this->m_Stencil->SetScalarComponentFromDouble + ( idx[ 0 ], idx[ 1 ], 0, 3, alpha ); + this->m_Stencil->Modified( ); + this->m_StencilActor->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image2DObserver< F, R >:: +Render( ) +{ + if( this->m_RenderWindow != NULL ) + this->m_RenderWindow->Render( ); +} + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image2DObserver< F, R >:: +Execute( const itk::Object* c, const itk::EventObject& e ) +{ + typedef typename TImage::IndexType TIndex; + typedef typename TImage::PointType TPoint; + typedef typename TFilter::TEvent TEvent; + typedef typename TFilter::TFrontEvent TFrontEvent; + typedef typename TFilter::TMarkEvent TMarkEvent; + typedef typename TFilter::TCollisionEvent TCollisionEvent; + typedef typename TFilter::TEndEvent TEndEvent; + + static unsigned char Colors[][4] = + { + { 0, 0, 127, 127 }, + { 0, 127, 127, 127 }, + { 127, 0, 127, 127 }, + { 127, 127, 0, 127 }, + { 0, 0, 63, 127 }, + { 0, 63, 63, 127 }, + { 63, 0, 63, 127 }, + { 63, 63, 0, 127 }, + { 63, 63, 127, 127 }, + { 63, 127, 127, 127 }, + { 127, 63, 127, 127 }, + { 127, 127, 63, 127 }, + { 127, 127, 63, 127 }, + { 127, 63, 63, 127 }, + { 63, 127, 63, 127 }, + { 63, 63, 127, 127 } + }; + + if( this->m_RenderWindow == NULL ) + return; + + const TEvent* evt = dynamic_cast< const TEvent* >( &e ); + TFrontEvent fevt; + TMarkEvent mevt; + TCollisionEvent cevt; + TEndEvent eevt; + if( evt != NULL ) + { + if( mevt.CheckEvent( evt ) ) + this->SetPixel( + evt->Node.Vertex, + Colors[ evt->Node.FrontId ][ 0 ], + Colors[ evt->Node.FrontId ][ 1 ], + Colors[ evt->Node.FrontId ][ 2 ], + Colors[ evt->Node.FrontId ][ 3 ] + ); + else if( fevt.CheckEvent( evt ) ) + this->SetPixel( evt->Node.Vertex, 255, 0, 0, 255 ); + if( cevt.CheckEvent( evt ) ) + this->SetPixel( evt->Node.Vertex, 100, 50, 25, 255 ); + + // Update visualization + if( this->m_Number == 0 || this->m_Percentage < 0 ) + this->Render( ); + this->m_Number++; + this->m_Number %= this->m_Percentage; + + if( eevt.CheckEvent( evt ) ) + { + if( this->m_RenderWindow != NULL ) + { + /* TODO + vtkRenderer* ren = + this->m_RenderWindow->GetRenderers( )->GetFirstRenderer( ); + ren->RemoveActor( this->m_StencilActor ); + */ + this->Render( ); + + } // fi + + } // fi + } + else + { + // TODO: Remove all! + this->Render( ); + + } // fi +} + +#endif // __FPA__VTK__IMAGE2DOBSERVER__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/Image3DObserver.h b/lib/fpa/VTK/Image3DObserver.h new file mode 100644 index 0000000..860946d --- /dev/null +++ b/lib/fpa/VTK/Image3DObserver.h @@ -0,0 +1,93 @@ +#ifndef __FPA__VTK__IMAGE3DOBSERVER__H__ +#define __FPA__VTK__IMAGE3DOBSERVER__H__ + +#include + +#include + +#include +#include +#include +#include + +namespace fpa +{ + namespace VTK + { + /** + */ + template< class F, class R > + class Image3DObserver + : public itk::Command + { + public: + typedef Image3DObserver Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef F TFilter; + typedef R TRenderWindow; + typedef typename TFilter::TInputImage TImage; + typedef typename TFilter::TVertex TVertex; + + typedef std::pair< unsigned long, unsigned long > TVertexIds; + + struct TVertexCmp + { + bool operator()( const TVertex& a, const TVertex& b ) const + { + unsigned int d = 0; + while( d < TVertex::Dimension && a[ d ] == b[ d ] ) + d++; + if( d < TVertex::Dimension ) + return( a[ d ] < b[ d ] ); + else + return( false ); + } + }; + typedef std::map< TVertex, TVertexIds, TVertexCmp > TVertices; + + public: + itkNewMacro( Self ); + itkTypeMacro( Image3DObserver, itkCommand ); + + public: + void SetImage( const TImage* img, R* rw ); + void Render( ); + void Execute( itk::Object* c, const itk::EventObject& e ) + { this->Execute( ( const itk::Object* )( c ), e ); } + void Execute( const itk::Object* c, const itk::EventObject& e ); + + protected: + Image3DObserver( ); + virtual ~Image3DObserver( ) + { } + + private: + Image3DObserver( const Self& ); // Not impl. + void operator=( const Self& ); // Not impl. + + protected: + typename TImage::ConstPointer m_Image; + + R* m_RenderWindow; + unsigned long m_Number; + unsigned long m_Percentage; + + vtkSmartPointer< vtkPolyData > m_PolyData; + vtkSmartPointer< vtkPolyDataMapper > m_PolyDataMapper; + vtkSmartPointer< vtkActor > m_PolyDataActor; + + TVertices m_Vertices; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__VTK__IMAGE3DOBSERVER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/Image3DObserver.hxx b/lib/fpa/VTK/Image3DObserver.hxx new file mode 100644 index 0000000..e79bbbe --- /dev/null +++ b/lib/fpa/VTK/Image3DObserver.hxx @@ -0,0 +1,203 @@ +#ifndef __FPA__VTK__IMAGE3DOBSERVER__HXX__ +#define __FPA__VTK__IMAGE3DOBSERVER__HXX__ + +#include +#include +#include +#include +#include +/* + #include + #include + #include + #include +*/ + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image3DObserver< F, R >:: +SetImage( const TImage* img, R* rw ) +{ + this->m_Image = img; + this->m_RenderWindow = rw; + this->m_Vertices.clear( ); + + /* + unsigned int minD = TImage::ImageDimension; + minD = ( minD < 3 )? minD: 3; + + int e[ 6 ] = { 0 }; + typename TImage::RegionType reg = this->m_Image->GetRequestedRegion( ); + for( unsigned int i = 0; i < minD; i++ ) + { + e[ ( i << 1 ) + 0 ] = reg.GetIndex( )[ i ]; + e[ ( i << 1 ) + 1 ] = reg.GetIndex( )[ i ] + reg.GetSize( )[ i ] - 1; + + } // rof + + typename TImage::SpacingType spac = this->m_Image->GetSpacing( ); + double s[ 3 ] = { 1, 1, 1 }; + for( unsigned int i = 0; i < minD; i++ ) + s[ i ] = double( spac[ i ] ); + + typename TImage::PointType orig = this->m_Image->GetOrigin( ); + double o[ 3 ] = { 0 }; + for( unsigned int i = 0; i < minD; i++ ) + o[ i ] = double( orig[ i ] ); + + this->m_Stencil->SetExtent( e ); + this->m_Stencil->SetSpacing( s ); + this->m_Stencil->SetOrigin( o ); + this->m_Stencil->AllocateScalars( VTK_UNSIGNED_CHAR, 4 ); + for( unsigned int i = 0; i < 3; i++ ) + this->m_Stencil->GetPointData( )-> + GetScalars( )->FillComponent( i, 255 ); + this->m_Stencil->GetPointData( )->GetScalars( )->FillComponent( 3, 0 ); + + this->m_StencilActor->SetInputData( this->m_Stencil ); + this->m_StencilActor->InterpolateOff( ); + */ + + double nPix = + double( this->m_Image->GetRequestedRegion( ).GetNumberOfPixels( ) ); + this->m_Percentage = ( unsigned long )( nPix * 0.001 ); + this->m_Number = 0; + + if( this->m_RenderWindow != NULL ) + { + vtkRenderer* ren = + this->m_RenderWindow->GetRenderers( )->GetFirstRenderer( ); + ren->AddActor( this->m_PolyDataActor ); + this->Render( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image3DObserver< F, R >:: +Render( ) +{ + if( this->m_RenderWindow != NULL ) + this->m_RenderWindow->Render( ); +} + +// ------------------------------------------------------------------------- +template< class F, class R > +void fpa::VTK::Image3DObserver< F, R >:: +Execute( const itk::Object* c, const itk::EventObject& e ) +{ + typedef typename TImage::IndexType TIndex; + typedef typename TImage::PointType TPoint; + typedef typename TFilter::TEvent TEvent; + typedef typename TFilter::TFrontEvent TFrontEvent; + typedef typename TFilter::TMarkEvent TMarkEvent; + typedef typename TFilter::TCollisionEvent TCollisionEvent; + typedef typename TFilter::TEndEvent TEndEvent; + + if( this->m_RenderWindow == NULL ) + return; + + const TEvent* evt = dynamic_cast< const TEvent* >( &e ); + TFrontEvent fevt; + TMarkEvent mevt; + TCollisionEvent cevt; + TEndEvent eevt; + if( evt != NULL ) + { + typename TImage::PointType pnt; + this->m_Image->TransformIndexToPhysicalPoint( evt->Node.Vertex, pnt ); + + if( fevt.CheckEvent( evt ) ) + { + if( + this->m_Vertices.find( evt->Node.Vertex ) == this->m_Vertices.end( ) + ) + { + unsigned long pId = this->m_PolyData->GetPoints( )->InsertNextPoint( + double( pnt[ 0 ] ), double( pnt[ 1 ] ), double( pnt[ 2 ] ) + ); + unsigned long cId = this->m_PolyData->GetVerts( )->InsertNextCell( 1 ); + this->m_PolyData->GetVerts( )->InsertCellPoint( pId ); + this->m_PolyData->Modified( ); + this->m_Vertices[ evt->Node.Vertex ] = TVertexIds( pId, cId ); + + } // rof + } + else if( mevt.CheckEvent( evt ) ) + { + /* + typename TVertices::iterator vIt = + this->m_Vertices.find( evt->Node.Vertex ); + if( vIt != this->m_Vertices.end( ) ) + { + std::cout << "Erase " << evt->Node.Vertex << " " << vIt->second.second << std::endl; + } // fi + */ + } // fi + + if( cevt.CheckEvent( evt ) ) + { + // std::cout << "Collision"; + } // fi + /* + this->SetPixel( + evt->Node.Vertex, + Colors[ evt->Node.FrontId ][ 0 ], + Colors[ evt->Node.FrontId ][ 1 ], + Colors[ evt->Node.FrontId ][ 2 ], + Colors[ evt->Node.FrontId ][ 3 ] + ); + else if( mevt.CheckEvent( evt ) ) + this->SetPixel( evt->Node.Vertex, 255, 0, 0, 255 ); + + if( cevt.CheckEvent( evt ) ) + this->SetPixel( evt->Node.Vertex, 100, 50, 25, 255 ); + */ + + // Update visualization + if( this->m_Number == 0 || this->m_Percentage < 0 ) + this->Render( ); + this->m_Number++; + this->m_Number %= this->m_Percentage; + + if( eevt.CheckEvent( evt ) ) + { + if( this->m_RenderWindow != NULL ) + { + vtkRenderer* ren = + this->m_RenderWindow->GetRenderers( )->GetFirstRenderer( ); + ren->RemoveActor( this->m_PolyDataActor ); + this->Render( ); + + } // fi + + } // fi + + } // fi +} + +// ------------------------------------------------------------------------- +template< class F, class R > +fpa::VTK::Image3DObserver< F, R >:: +Image3DObserver( ) + : Superclass( ), + m_RenderWindow( NULL ), + m_Number( 0 ), + m_Percentage( 0 ) +{ + this->m_PolyData = vtkSmartPointer< vtkPolyData >::New( ); + this->m_PolyDataMapper = vtkSmartPointer< vtkPolyDataMapper >::New( ); + this->m_PolyDataActor = vtkSmartPointer< vtkActor >::New( ); + + this->m_PolyData->SetPoints( vtkPoints::New( ) ); + this->m_PolyData->SetVerts( vtkCellArray::New( ) ); + + this->m_PolyDataMapper->SetInputData( this->m_PolyData ); + this->m_PolyDataActor->SetMapper( this->m_PolyDataMapper ); + this->m_PolyDataActor->GetProperty( )->SetColor( 0, 0, 1 ); +} + +#endif // __FPA__VTK__IMAGE3DOBSERVER__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/ImageMPR.cxx b/lib/fpa/VTK/ImageMPR.cxx new file mode 100644 index 0000000..d435bb8 --- /dev/null +++ b/lib/fpa/VTK/ImageMPR.cxx @@ -0,0 +1,159 @@ +#include + +#include + +// ------------------------------------------------------------------------- +fpa::VTK::ImageMPR:: +ImageMPR( ) +{ + // vtkSmartPointer< vtkImageData > m_Image; + + this->m_Outline = vtkSmartPointer< vtkOutlineSource >::New( ); + this->m_OutlineMapper = vtkSmartPointer< vtkPolyDataMapper >::New( ); + this->m_OutlineActor = vtkSmartPointer< vtkActor >::New( ); + this->m_Picker = vtkSmartPointer< vtkCellPicker >::New( ); + this->m_WidgetX = vtkSmartPointer< vtkImagePlaneWidget >::New( ); + this->m_WidgetY = vtkSmartPointer< vtkImagePlaneWidget >::New( ); + this->m_WidgetZ = vtkSmartPointer< vtkImagePlaneWidget >::New( ); + this->m_Renderer = vtkSmartPointer< vtkRenderer >::New( ); + this->m_Window = vtkSmartPointer< vtkRenderWindow >::New( ); + this->m_Interactor = vtkSmartPointer< vtkRenderWindowInteractor >::New( ); +} + +// ------------------------------------------------------------------------- +fpa::VTK::ImageMPR:: +~ImageMPR( ) +{ +} + +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +SetImage( vtkImageData* image ) +{ + this->m_Image = image; + + // Outline + this->m_Outline->SetBounds( this->m_Image->GetBounds( ) ); + this->m_OutlineMapper-> + SetInputConnection( this->m_Outline->GetOutputPort( ) ); + this->m_OutlineActor->SetMapper( this->m_OutlineMapper ); + this->m_OutlineActor->GetProperty( )->SetColor( 1, 1, 1 ); + + // Local picker + this->m_Picker->SetTolerance( 0.005 ); + + // Image planes + this->m_WidgetX->DisplayTextOn( ); + this->m_WidgetX->SetInputData( this->m_Image ); + this->m_WidgetX->SetPlaneOrientationToXAxes( ); + this->m_WidgetX->SetSliceIndex( this->m_Image->GetExtent( )[ 0 ] ); + this->m_WidgetX->SetPicker( this->m_Picker ); + this->m_WidgetX->SetKeyPressActivationValue( 'x' ); + this->m_WidgetX->GetPlaneProperty( )->SetLineWidth( 3 ); + this->m_WidgetX->GetPlaneProperty( )->SetColor( 0, 1, 1 ); + + this->m_WidgetY->DisplayTextOn( ); + this->m_WidgetY->SetInputData( this->m_Image ); + this->m_WidgetY->SetPlaneOrientationToYAxes( ); + this->m_WidgetY->SetSliceIndex( this->m_Image->GetExtent( )[ 2 ] ); + this->m_WidgetY->SetPicker( this->m_Picker ); + this->m_WidgetY->SetKeyPressActivationValue( 'y' ); + this->m_WidgetY->GetPlaneProperty( )->SetColor( 1, 0, 1 ); + this->m_WidgetY->GetPlaneProperty( )->SetLineWidth( 3 ); + this->m_WidgetY->SetLookupTable( this->m_WidgetX->GetLookupTable( ) ); + + this->m_WidgetZ->DisplayTextOn( ); + this->m_WidgetZ->SetInputData( this->m_Image ); + this->m_WidgetZ->SetPlaneOrientationToZAxes( ); + this->m_WidgetZ->SetSliceIndex( this->m_Image->GetExtent( )[ 4 ] ); + this->m_WidgetZ->SetPicker( this->m_Picker ); + this->m_WidgetZ->SetKeyPressActivationValue( 'z' ); + this->m_WidgetZ->GetPlaneProperty( )->SetColor( 1, 1, 0 ); + this->m_WidgetZ->GetPlaneProperty( )->SetLineWidth( 3 ); + this->m_WidgetZ->SetLookupTable( this->m_WidgetX->GetLookupTable( ) ); + + // Rendering stuff + this->m_Window->AddRenderer( this->m_Renderer ); + this->m_Interactor->SetRenderWindow( this->m_Window ); + + // Command + double spac[ 3 ]; + this->m_Image->GetSpacing( spac ); + double min_spacing = spac[ 0 ]; + for( unsigned int d = 1; d < 3; d++ ) + min_spacing = ( spac[ d ] < min_spacing )? spac[ d ]: min_spacing; + + vtkInteractorStyleSwitch* iswitch = + dynamic_cast< vtkInteractorStyleSwitch* >( + this->m_Interactor->GetInteractorStyle( ) + ); + if( iswitch != NULL ) + iswitch->SetCurrentStyleToTrackballCamera( ); + + // Add actors + this->m_Renderer->AddActor( this->m_OutlineActor ); + + // Prepare widgets + this->m_WidgetX->SetInteractor( this->m_Interactor ); + this->m_WidgetY->SetInteractor( this->m_Interactor ); + this->m_WidgetZ->SetInteractor( this->m_Interactor ); +} + +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +SetBackground( double r, double g, double b ) +{ + this->m_Renderer->SetBackground( r, g, b ); +} + +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +SetSize( unsigned int w, unsigned int h ) +{ + this->m_Window->SetSize( w, h ); +} + +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +AddPolyData( vtkPolyData* pd, double r, double g, double b ) +{ + unsigned int i = this->m_PolyDatas.size( ); + + this->m_PolyDatas.push_back( pd ); + this->m_Mappers.push_back( vtkSmartPointer< vtkPolyDataMapper >::New( ) ); + this->m_Actors.push_back( vtkSmartPointer< vtkActor >::New( ) ); + + this->m_Mappers[ i ]->SetInputData( pd ); + this->m_Actors[ i ]->SetMapper( this->m_Mappers[ i ] ); + this->m_Actors[ i ]->GetProperty( )->SetColor( r, g, b ); + this->m_Renderer->AddActor( this->m_Actors[ i ] ); +} + +// ------------------------------------------------------------------------- +vtkRenderWindow* fpa::VTK::ImageMPR:: +GetWindow( ) const +{ + return( this->m_Window ); +} + +// ------------------------------------------------------------------------- +vtkRenderer* fpa::VTK::ImageMPR:: +GetRenderer( ) const +{ + return( this->m_Renderer ); +} + +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +Start( ) +{ + this->m_WidgetX->On( ); + this->m_WidgetY->On( ); + this->m_WidgetZ->On( ); + + this->m_Renderer->ResetCamera( ); + this->m_Interactor->Initialize( ); + this->m_Interactor->Start( ); +} + +// eof - $RCSfile$ diff --git a/lib/fpa/VTK/ImageMPR.h b/lib/fpa/VTK/ImageMPR.h new file mode 100644 index 0000000..e7a65d3 --- /dev/null +++ b/lib/fpa/VTK/ImageMPR.h @@ -0,0 +1,67 @@ +#ifndef __FPA__VTK__IMAGEMPR__H__ +#define __FPA__VTK__IMAGEMPR__H__ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace VTK + { + /** + */ + class FrontAlgorithms_EXPORT ImageMPR + { + public: + ImageMPR( ); + virtual ~ImageMPR( ); + + void SetImage( vtkImageData* image ); + void SetBackground( double r, double g, double b ); + void SetSize( unsigned int w, unsigned int h ); + + void AddPolyData( vtkPolyData* pd, double r, double g, double b ); + + vtkRenderWindow* GetWindow( ) const; + vtkRenderer* GetRenderer( ) const; + + void Start( ); + + protected: + vtkSmartPointer< vtkImageData > m_Image; + vtkSmartPointer< vtkOutlineSource > m_Outline; + vtkSmartPointer< vtkPolyDataMapper > m_OutlineMapper; + vtkSmartPointer< vtkActor > m_OutlineActor; + vtkSmartPointer< vtkCellPicker > m_Picker; + vtkSmartPointer< vtkImagePlaneWidget > m_WidgetX; + vtkSmartPointer< vtkImagePlaneWidget > m_WidgetY; + vtkSmartPointer< vtkImagePlaneWidget > m_WidgetZ; + vtkSmartPointer< vtkRenderer > m_Renderer; + vtkSmartPointer< vtkRenderWindow > m_Window; + vtkSmartPointer< vtkRenderWindowInteractor > m_Interactor; + + std::vector< vtkSmartPointer< vtkPolyData > > m_PolyDatas; + std::vector< vtkSmartPointer< vtkPolyDataMapper > > m_Mappers; + std::vector< vtkSmartPointer< vtkActor > > m_Actors; + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__VTK__IMAGEMPR__H__ + +// eof - $RCSfile$ -- 2.45.1