#include <chrono>
#include <map>
+
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
-#include <fpa/Filters/Image/Mori.h>
+#include <itkRegionOfInterestImageFilter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+
+#include <fpa/Filters/BaseMarksInterface.h>
+#include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
+#include <fpa/Filters/Image/DefaultTraits.h>
+#include <fpa/Filters/Image/RegionGrow.h>
+#include <fpa/Functors/RegionGrow/BinaryThreshold.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;
+
+/**
+ */
+class MoriLabellingTraits
+ : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
+{
+public:
+ typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass;
+ typedef typename Superclass::TInternalTraits TInternalTraits;
+ typedef typename Superclass::TMarksImage TMarksImage;
+ typedef typename Superclass::TFilterInterface TFilterInterface;
+
+ typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
+ typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
+};
+
+/**
+ */
+class MoriLabelling
+ : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
+{
+public:
+ typedef MoriLabellingTraits TTraits;
+ fpaTraitsMacro( TTraits );
+
+ typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass;
+ typedef MoriLabelling Self;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+ typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor;
+
+public:
+ itkNewMacro( Self );
+ itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
+
+ itkGetConstMacro( LastThreshold, TInputValue );
+ itkSetMacro( LastThreshold, TInputValue );
+
+ itkGetConstMacro( MinVertex, TVertex );
+ itkGetConstMacro( MaxVertex, TVertex );
+
+ fpaFilterInputMacro( InputLabels, TLabels );
+ fpaFilterInputMacro( InputVesselness, TScalarImage );
+
+public:
+ TInputValue GetUpperThreshold( ) const
+ {
+ return( this->m_Functor->GetUpperThreshold( ) );
+ }
+ void SetUpperThreshold( TInputValue t )
+ {
+ this->m_Functor->SetUpperThreshold( t );
+ }
+
+protected:
+ MoriLabelling( )
+ : Superclass( ),
+ 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 );
+ }
+
+ virtual ~MoriLabelling( )
+ {
+ }
+
+ virtual const itk::DataObject* _GetReferenceInput( ) const override
+ {
+ return( this->GetInputLabels( ) );
+ }
+ virtual void _BeforeGenerateData( ) override
+ {
+ 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
+ {
+ this->Superclass::_PostComputeOutputValue( n );
+ if( n.Value == this->GetInsideValue( ) )
+ {
+ 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 )
+ {
+ if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
+ n.Value = this->GetInsideValue( );
+ else
+ n.Value = 0;
+
+ } // fi
+
+ if( !( this->m_FirstVertex ) )
+ {
+ for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
+ {
+ if( n.Vertex[ d ] < this->m_MinVertex[ d ] )
+ this->m_MinVertex[ d ] = n.Vertex[ d ];
+ if( this->m_MaxVertex[ d ] < n.Vertex[ d ] )
+ this->m_MaxVertex[ d ] = n.Vertex[ d ];
+
+ } // rof
+ }
+ else
+ {
+ this->m_MinVertex = n.Vertex;
+ this->m_MaxVertex = n.Vertex;
+ this->m_FirstVertex = false;
+
+ } // fi
+
+ } // fi
+ }
+
+private:
+ // Purposely not implemented.
+ MoriLabelling( const Self& other );
+ Self& operator=( const Self& other );
+
+protected:
+ TFunctor::Pointer m_Functor;
+ TInputValue m_LastThreshold;
+ bool m_FirstVertex;
+ TVertex m_MinVertex;
+ TVertex m_MaxVertex;
+
+ TScalar m_MaxVesselness;
+ TScalar m_VesselnessThr;
+};
// -------------------------------------------------------------------------
double MeasureTime( itk::ProcessObject* f )
// -------------------------------------------------------------------------
template< class _TImage >
-void WriteImage( const _TImage*image, const std::string& fname )
+void WriteImage( const _TImage* image, const std::string& fname )
{
typedef itk::ImageFileWriter< _TImage > _TWriter;
typename _TWriter::Pointer writer = _TWriter::New( );
std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
}
+// -------------------------------------------------------------------------
+template< class _TImagePtr, class _TRegion >
+void ROI(
+ _TImagePtr& output,
+ const typename _TImagePtr::ObjectType* input,
+ const _TRegion& roi
+ )
+{
+ typedef typename _TImagePtr::ObjectType _TImage;
+ typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
+
+ typename _TFilter::Pointer filter = _TFilter::New( );
+ filter->SetInput( input );
+ filter->SetRegionOfInterest( roi );
+ double t = MeasureTime( filter );
+ std::cout << "ROI computed in " << t << " s" << std::endl;
+ output = filter->GetOutput( );
+ output->DisconnectPipeline( );
+}
+
// -------------------------------------------------------------------------
bool ParseArgs(
std::map< std::string, std::string >& args, int argc, char* argv[]
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" ] );
- // Mori segmentation
+ // 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( ) ) )
+ );
/* TODO
- typedef fpa::Filters::Image::Mori< TImage, TLabels > TMori;
- TMori::Pointer mori = TMori::New( );
- mori->SetInput( input_image );
- mori->SetSeed( seed );
- mori->SetInsideValue( 1 );
- mori->SetOutsideValue( 0 );
- mori->SetMinimumThreshold(
- TPixel( std::atof( args[ "minimum_threshold" ].c_str( ) ) )
- );
- mori->SetSignalKernelSize(
- std::atoi( args[ "signal_kernel_size" ].c_str( ) )
- );
- mori->SetSignalThreshold(
- std::atof( args[ "signal_threshold" ].c_str( ) )
- );
- mori->SetSignalInfluence(
- std::atof( args[ "signal_influence" ].c_str( ) )
- );
- mori->SetThresholds(
- TPixel( std::atof( args[ "lower_threshold" ].c_str( ) ) ),
- TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ),
- TPixel( std::atof( args[ "delta_threshold" ].c_str( ) ) )
- );
- double t = MeasureTime( mori );
- std::cout << "Mori executed in " << t << " s" << std::endl;
- WriteImage( mori->GetOutput( ), args[ "out" ] );
-
- std::map< std::string, std::string >::const_iterator i =
- args.find( "out_signal" );
- if( i != args.end( ) )
- {
- std::stringstream signal;
- unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
- signal << "# nThr = " << nthr << std::endl;
- signal << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl;
- for( unsigned long j = 0; j < nthr; ++j )
- {
- typename TMori::TPeak p;
- double x, y;
- mori->GetSignalValues( j, x, y, p );
- signal << x << " " << y << std::endl;
-
- } // rof
- std::ofstream signals_str( i->second.c_str( ) );
- signals_str << signal.str( );
- signals_str.close( );
-
- } // fi
+ labelling->SetLastThreshold( last_thr );
*/
+ double t = MeasureTime( labelling );
+ std::cout << "Labelling executed in " << t << " s" << std::endl;
+
+ std::map< std::string, std::string >::const_iterator i =
+ args.find( "roi" );
+ if( i != args.end( ) )
+ {
+ TImage::IndexType minV = labelling->GetMinVertex( );
+ TImage::IndexType maxV = labelling->GetMaxVertex( );
+ TImage::SizeType roiS;
+ for( unsigned d = 0; d < TImage::ImageDimension; ++d )
+ roiS[ d ] = maxV[ d ] - minV[ d ] + 1;
+ TImage::RegionType roi;
+ roi.SetIndex( minV );
+ roi.SetSize( roiS );
+
+ // ROI input image
+ TImage::Pointer input_image_roi;
+ ROI( input_image_roi, input_image.GetPointer( ), roi );
+ WriteImage( input_image_roi.GetPointer( ), i->second );
+
+ // ROI output image
+ TLabels::Pointer output_roi;
+ ROI( output_roi, labelling->GetOutput( ), roi );
+ WriteImage( output_roi.GetPointer( ), args[ "out" ] );
+ }
+ else
+ WriteImage( labelling->GetOutput( ), args[ "out" ] );
}
else
return( 1 );
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff