${PROJECT_SOURCE_DIR}/lib
${PROJECT_BINARY_DIR}/lib
)
-subdirs(CTArteries)
+subdirs(CTArteries CTBronchi)
## eof - $RCSfile$
+## =========================================================================
+## @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+## =========================================================================
+
option(fpa_BUILD_CTBronchi "Build bronchi analysis from CT images applications?" OFF)
if(fpa_BUILD_CTBronchi)
- configure_file(
- CTBronchi_process.sh
- ${PROJECT_BINARY_DIR}/CTBronchi_process.sh
- COPYONLY
+ set(_pfx fpa_CTBronchi_)
+ set(
+ _examples
+ Vesselness
+ MoriSegmentation
+ MoriLabelling
)
- include_directories(
- ${PROJECT_SOURCE_DIR}/lib
- ${PROJECT_BINARY_DIR}/lib
- ${PROJECT_SOURCE_DIR}/appli
- ${PROJECT_BINARY_DIR}/appli
- )
-set(_pfx fpa_CTBronchi_)
-set(
- _examples
- MoriSegmentation
- MoriLabelling
- )
-foreach(_e ${_examples})
- add_executable(${_pfx}${_e} ${_e}.cxx)
- target_link_libraries(${_pfx}${_e} fpa)
-endforeach(_e)
+ foreach(_e ${_examples})
+ BuildApplication(
+ ${_pfx}${_e}
+ SOURCE ${_e}.cxx
+ INSTALL
+ RECURRENT
+ LINKS fpa
+ )
+ endforeach(_e)
endif(fpa_BUILD_CTBronchi)
## eof - $RCSfile$
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkRegionOfInterestImageFilter.h>
+#include <itkMinimumMaximumImageCalculator.h>
#include <fpa/Filters/BaseMarksInterface.h>
#include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
-typedef short TPixel;
-typedef unsigned char TLabel;
-typedef itk::Image< TPixel, Dim > TImage;
-typedef itk::Image< TLabel, Dim > TLabels;
+typedef short TPixel;
+typedef unsigned char TLabel;
+typedef itk::NumericTraits< TPixel >::RealType TScalar;
+typedef itk::Image< TPixel, Dim > TImage;
+typedef itk::Image< TLabel, Dim > TLabels;
+typedef itk::Image< TScalar, Dim > TScalarImage;
/**
*/
itkGetConstMacro( MaxVertex, TVertex );
fpaFilterInputMacro( InputLabels, TLabels );
+ fpaFilterInputMacro( InputVesselness, TScalarImage );
public:
TInputValue GetUpperThreshold( ) const
protected:
MoriLabelling( )
: Superclass( ),
- m_LastThreshold( TInputValue( 0 ) )
+ m_LastThreshold( TInputValue( 0 ) ),
+ m_VesselnessThr( TScalar( 0.05 ) )
{
fpaFilterInputConfigureMacro( InputLabels, TLabels );
+ fpaFilterInputConfigureMacro( InputVesselness, TScalarImage );
this->m_Functor = TFunctor::New( );
this->SetPredicate( this->m_Functor );
}
{
this->Superclass::_BeforeGenerateData( );
this->m_FirstVertex = true;
+
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ _TMinMax::Pointer minMax = _TMinMax::New( );
+ minMax->SetImage( this->GetInputVesselness( ) );
+ minMax->Compute( );
+ this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( );
}
virtual void _PostComputeOutputValue( TNode& n ) override
{
const TImage* input = this->GetInput( );
const TLabels* labels = this->GetInputLabels( );
+ const TScalarImage* vesselness = this->GetInputVesselness( );
double x = input->GetPixel( n.Vertex );
/* TODO
double a = std::fabs( x - double( this->m_Functor->GetUpperThreshold( ) ) );
double b = std::fabs( x - double( this->m_LastThreshold ) );
*/
- if( labels->GetPixel( n.Vertex ) == 0 /* && b < a*/ )
- n.Value = 0;
+ if( labels->GetPixel( n.Vertex ) == 0 )
+ {
+ if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
+ n.Value = this->GetInsideValue( );
+ else
+ n.Value = 0;
+
+ } // fi
if( !( this->m_FirstVertex ) )
{
bool m_FirstVertex;
TVertex m_MinVertex;
TVertex m_MaxVertex;
+
+ TScalar m_MaxVesselness;
+ TScalar m_VesselnessThr;
};
// -------------------------------------------------------------------------
mandatory.insert( "in" );
mandatory.insert( "out" );
mandatory.insert( "labels" );
+ mandatory.insert( "vesselness" );
- args[ "upper_threshold" ] = "-600";
+ args[ "upper_threshold" ] = "-650";
int i = 1;
while( i < argc )
<< "\t-in filename" << std::endl
<< "\t-out filename" << std::endl
<< "\t-labels filename" << std::endl
+ << "\t-vesselness filename" << std::endl
<< "\t[-roi] filename" << std::endl
<< "\t[-upper_threshold value]" << std::endl;
return( false );
TLabels::Pointer input_labels;
ReadImage( input_labels, args[ "labels" ] );
+ // Read vesselness image
+ TScalarImage::Pointer input_vesselness;
+ ReadImage( input_vesselness, args[ "vesselness" ] );
+
// Mori labelling
MoriLabelling::Pointer labelling = MoriLabelling::New( );
labelling->SetInput( input_image );
labelling->SetInputLabels( input_labels );
+ labelling->SetInputVesselness( input_vesselness );
labelling->SetOutsideValue( 2 );
+ labelling->SetFillValue( 2 );
labelling->SetInsideValue( 1 );
labelling->SetUpperThreshold(
TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) )
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#include <chrono>
+#include <map>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkInvertIntensityImageFilter.h>
+#include <itkHessianRecursiveGaussianImageFilter.h>
+#include <itkHessian3DToVesselnessMeasureImageFilter.h>
+
+// -------------------------------------------------------------------------
+const unsigned int Dim = 3;
+typedef short TPixel;
+typedef itk::NumericTraits< TPixel >::RealType TScalar;
+typedef itk::Image< TPixel, Dim > TImage;
+typedef itk::Image< TScalar, Dim > TScalarImage;
+
+// -------------------------------------------------------------------------
+double MeasureTime( itk::ProcessObject* f )
+{
+ std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
+ std::chrono::duration< double > t;
+ s = std::chrono::high_resolution_clock::now( );
+ f->Update( );
+ e = std::chrono::high_resolution_clock::now( );
+ t = e - s;
+ return( t.count( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImagePtr >
+void ReadImage( _TImagePtr& image, const std::string& fname )
+{
+ typedef typename _TImagePtr::ObjectType _TImage;
+ typedef itk::ImageFileReader< _TImage > _TReader;
+ typename _TReader::Pointer reader = _TReader::New( );
+ reader->SetFileName( fname );
+ double t = MeasureTime( reader );
+ std::cout << "Read " << fname << " in " << t << " s" << std::endl;
+ image = reader->GetOutput( );
+ image->DisconnectPipeline( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void WriteImage( const _TImage* image, const std::string& fname )
+{
+ typedef itk::ImageFileWriter< _TImage > _TWriter;
+ typename _TWriter::Pointer writer = _TWriter::New( );
+ writer->SetFileName( fname );
+ writer->SetInput( image );
+ double t = MeasureTime( writer );
+ std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
+}
+
+// -------------------------------------------------------------------------
+bool ParseArgs(
+ std::map< std::string, std::string >& args, int argc, char* argv[]
+ )
+{
+ std::set< std::string > mandatory;
+ mandatory.insert( "in" );
+ mandatory.insert( "out" );
+
+ args[ "sigma" ] = "0.5";
+ args[ "alpha1" ] = "0.5";
+ args[ "alpha2" ] = "2";
+
+ int i = 1;
+ while( i < argc )
+ {
+ std::string cmd = argv[ i ] + 1;
+ args[ cmd ] = argv[ i + 1 ];
+ i += 2;
+
+ } // elihw
+
+ bool complete = true;
+ for( std::string t: mandatory )
+ complete &= ( args.find( t ) != args.end( ) );
+
+ if( !complete )
+ {
+ std::cerr
+ << "Usage: " << argv[ 0 ] << std::endl
+ << "\t-in filename" << std::endl
+ << "\t-out filename" << std::endl
+ << "\t[-sigma val]" << std::endl
+ << "\t[-alpha1 val]" << std::endl
+ << "\t[-alpha2 val]" << std::endl;
+ return( false );
+
+ } // fi
+ return( true );
+}
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+ std::map< std::string, std::string > args;
+ try
+ {
+ if( ParseArgs( args, argc, argv ) )
+ {
+ // Read input image
+ TImage::Pointer input_image;
+ ReadImage( input_image, args[ "in" ] );
+
+ // Vesselness
+ typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax;
+ typedef itk::InvertIntensityImageFilter< TImage > TInverter;
+ typedef itk::HessianRecursiveGaussianImageFilter< TImage > THessian;
+ typedef itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > TVesselness;
+
+ TMinMax::Pointer minMax = TMinMax::New( );
+ minMax->SetImage( input_image );
+ minMax->Compute( );
+
+ TInverter::Pointer inverter = TInverter::New( );
+ inverter->SetInput( input_image );
+ inverter->SetMaximum( minMax->GetMaximum( ) );
+ double t = MeasureTime( inverter );
+ std::cout << "Inversion executed in " << t << " s" << std::endl;
+
+ THessian::Pointer hessian = THessian::New( );
+ hessian->SetInput( inverter->GetOutput( ) );
+ hessian->SetSigma(
+ std::atof( args[ "sigma" ].c_str( ) )
+ );
+ t = MeasureTime( hessian );
+ std::cout << "Hessian executed in " << t << " s" << std::endl;
+
+ TVesselness::Pointer vesselness = TVesselness::New( );
+ vesselness->SetInput( hessian->GetOutput( ) );
+ vesselness->SetAlpha1(
+ std::atof( args[ "alpha1" ].c_str( ) )
+ );
+ vesselness->SetAlpha2(
+ std::atof( args[ "alpha2" ].c_str( ) )
+ );
+ t = MeasureTime( vesselness );
+ std::cout << "Vesselness executed in " << t << " s" << std::endl;
+
+ WriteImage( vesselness->GetOutput( ), args[ "out" ] );
+ }
+ else
+ return( 1 );
+ }
+ catch( std::exception& err )
+ {
+ std::cerr
+ << "===============================" << std::endl
+ << "Error caught: " << std::endl
+ << err.what( )
+ << "===============================" << std::endl
+ << std::endl;
+ return( 1 );
+
+ } // yrt
+ return( 0 );
+}
+
+// eof - $RCSfile$
#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
// -------------------------------------------------------------------------
-const unsigned int Dim = 2;
-typedef unsigned char TPixel;
+const unsigned int Dim = 3;
+typedef short TPixel;
typedef unsigned char TLabel;
-typedef float TCost;
+typedef double TCost;
typedef itk::Image< TPixel, Dim > TImage;
typedef itk::Image< TLabel, Dim > TLabelImage;
TValue d = TValue( image->GetPixel( v ) );
d -= TValue( image->GetPixel( p ) );
d /= this->m_Beta;
- d *= d;
- if( this->m_TreatAsWeight ) d = std::exp( d ) - TValue( 1 );
- else d = std::exp( -d );
+ if( this->m_TreatAsWeight ) d = std::exp( d * d ) - TValue( 1 );
+ else d = std::exp( -std::fabs( d ) );
if( d < this->m_Epsilon ) return( this->m_Epsilon );
else return( d );
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff