#include <string>
#include <tclap/CmdLine.h>
#include <itkImage.h>
-#include "MoriLabelling.h"
+#include <fpa/Filters/Image/RandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
#include "Functions.h"
// -------------------------------------------------------------------------
typedef itk::NumericTraits< TPixel >::RealType TScalar;
typedef itk::Image< TPixel, Dim > TImage;
typedef itk::Image< TLabel, Dim > TLabels;
-typedef itk::Image< TScalar, Dim > TScalarImage;
// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
typedef TCLAP::ValueArg< std::string > _TStringArg;
- typedef TCLAP::ValueArg< double > _TRealArg;
- typedef TCLAP::ValueArg< TPixel > _TPixelArg;
+ typedef TCLAP::ValueArg< TScalar > _TRealArg;
// Parse input line
_TStringArg in( "i", "input", "Input image", true, "", "file" );
_TStringArg labels( "l", "labels", "Input labels", true, "", "file" );
- _TStringArg vesselness( "v", "vesselness", "Input vesselness", true, "", "file" );
_TStringArg out( "o", "output", "Output image", true, "", "file" );
- _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold", false, 0.05, "value (0.05)" );
- _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -400, "value (-400)" );
+ _TRealArg beta( "b", "beta", "Beta", false, 2.5, "value (2.5)" );
+ _TRealArg eps( "e", "epsilon", "Epsilon", false, 1e-5, "value (1e-5)" );
try
{
- TCLAP::CmdLine cmd( "Labelling", ' ', "1.0.0" );
- cmd.add( uThr );
- cmd.add( vThr );
+ TCLAP::CmdLine cmd( "FastRandomWalker", ' ', "1.0.0" );
+ cmd.add( eps );
+ cmd.add( beta );
cmd.add( out );
- cmd.add( vesselness );
cmd.add( labels );
cmd.add( in );
cmd.parse( argc, argv );
TLabels::Pointer input_labels;
CTBronchi::ReadImage( input_labels, labels.getValue( ) );
- // Read input vesselness
- TScalarImage::Pointer input_vesselness;
- CTBronchi::ReadImage( input_vesselness, vesselness.getValue( ) );
+ // Prepare weight functor
+ typedef fpa::Functors::Dijkstra::Image::Gaussian< TImage, TScalar > TWeight;
+ TWeight::Pointer weight = TWeight::New( );
+ weight->SetBeta( beta.getValue( ) );
+ weight->SetEpsilon( eps.getValue( ) );
- // Create labels
- typedef CTBronchi::MoriLabelling< TImage, TLabels, TScalarImage > _TLabelling;
- _TLabelling::Pointer labelling = _TLabelling::New( );
- labelling->SetInput( input_image );
- labelling->SetInputLabels( input_labels );
- labelling->SetInputVesselness( input_vesselness );
- labelling->SetVesselnessThreshold( vThr.getValue( ) );
- labelling->SetUpperThreshold( uThr.getValue( ) );
- labelling->SetInsideValue( 1 );
- labelling->SetOutsideValue( 2 );
- labelling->SetFillValue( 2 );
- double t = CTBronchi::MeasureTime( labelling );
- std::cout << "Labelling executed in " << t << " s" << std::endl;
+ // Random walk
+ typedef fpa::Filters::Image::RandomWalker< TImage, TLabels, TScalar > TFilter;
+ TFilter::Pointer filter = TFilter::New( );
+ filter->SetInputImage( input_image );
+ filter->SetInputLabels( input_labels );
+ filter->SetWeightFunction( weight );
+ double t = CTBronchi::MeasureTime( filter );
+ std::cout << "FastRandomWalk executed in " << t << " s" << std::endl;
// Write result
- CTBronchi::WriteImage( labelling->GetOutput( ), out.getValue( ) );
+ CTBronchi::WriteImage( filter->GetOutputLabels( ), out.getValue( ) );
}
catch( std::exception& err )
{
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __fpa__Common__SliceBySliceRandomWalker__h__
+#define __fpa__Common__SliceBySliceRandomWalker__h__
+
+#include <fpa/Config.h>
+#include <itkImageToImageFilter.h>
+
+namespace fpa
+{
+ namespace Common
+ {
+ /**
+ */
+ template< class _TImage, class _TLabels, class _TScalarImage >
+ class SliceBySliceRandomWalker
+ : public itk::ImageToImageFilter< _TImage, _TLabels >
+ {
+ public:
+ typedef _TImage TImage;
+ typedef _TLabels TLabels;
+ typedef _TScalarImage TScalarImage;
+
+ typedef SliceBySliceRandomWalker Self;
+ typedef itk::ImageToImageFilter< TImage, TLabels > Superclass;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ typedef typename TImage::PixelType TPixel;
+ typedef typename TLabels::PixelType TLabel;
+ typedef typename TScalarImage::PixelType TScalar;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro(
+ fpa::Common::SliceBySliceRandomWalker, itk::ImageToImageFilter
+ );
+
+ ivqITKInputMacro( InputLabels, TLabels );
+ ivqITKInputMacro( InputVesselness, TScalarImage );
+
+ protected:
+ SliceBySliceRandomWalker( );
+ virtual ~SliceBySliceRandomWalker( );
+
+ virtual void GenerateData( ) override;
+
+ void _Composite(
+ typename TScalarImage::Pointer& composite,
+ const TLabels* labels,
+ const TScalarImage* vesselness
+ );
+
+ template< class _TSlicePtr, class _TInput >
+ void _Slice(
+ _TSlicePtr& slice,
+ const _TInput* input,
+ typename _TInput::RegionType region
+ );
+
+ template< class _TBinaryTree, class _TData3D, class _TFusion >
+ void _GoDown(
+ _TBinaryTree& binaryTree,
+ const _TData3D& data3D,
+ const _TFusion& fusion,
+ int numSlices
+ );
+
+ template< class _TBinaryTree, class _TData3D, class _TFusion >
+ void _GoUp(
+ _TBinaryTree& binaryTree,
+ const _TData3D& data3D,
+ const _TFusion& fusion,
+ int numSlices
+ );
+
+ private:
+ // Purposely not implemented
+ SliceBySliceRandomWalker( const Self& other );
+ Self& operator=( const Self& other );
+ };
+
+ } // ecapseman
+
+} // ecapseman
+
+#ifndef ITK_MANUAL_INSTANTIATION
+# include <fpa/Common/SliceBySliceRandomWalker.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+#endif // __fpa__Common__SliceBySliceRandomWalker__h__
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__
+#define __fpa__Common__SliceBySliceRandomWalker__hxx__
+
+#include <vector>
+#include <itkExtractImageFilter.h>
+#include <itkImageRegionConstIterator.h>
+#include <itkImageRegionIterator.h>
+#include <itkJoinSeriesImageFilter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+SliceBySliceRandomWalker( )
+ : Superclass( )
+{
+ ivqITKInputConfigureMacro( InputLabels, TLabels );
+ ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+~SliceBySliceRandomWalker( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+void
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+GenerateData( )
+{
+ typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
+ typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
+ typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;
+
+ // Preare I/O objects
+ const TImage* input = this->GetInput( );
+ const TLabels* labels = this->GetInputLabels( );
+ const TScalarImage* vesselness = this->GetInputVesselness( );
+ this->AllocateOutputs( );
+
+ // Some values
+ typename TImage::RegionType region = input->GetRequestedRegion( );
+ typename TImage::SizeType regionIndex = region.GetIndex( );
+ typename TImage::SizeType sliceSize = region.GetSize( );
+ int numSlices = sliceSize[ 2 ];
+ sliceSize[ 2 ] = 0;
+
+ // Composite image
+ typename TScalarImage::Pointer composite = TScalarImage::New( );
+ composite->CopyInformation( vesselness );
+ composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
+ composite->Allocate( );
+ this->_Composite( composite, labels, vesselness );
+
+ // Extract slices
+ std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
+ std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
+ std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
+ for( int i = 0; i < numSlices; ++i )
+ {
+ // Prepare extraction region
+ typename TImage::IndexType sliceIndex = regionIndex;
+ sliceIndex[ 2 ] += i;
+ typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );
+
+ // Extract regions
+ typename _TSliceImage::Pointer input_slice;
+ typename _TSliceLabels::Pointer labels_slice;
+ typename _TSliceScalarImage::Pointer composite_slice;
+ this->_Slice( input_slice, input, sliceRegion );
+ this->_Slice( labels_slice, labels, sliceRegion );
+ this->_Slice( composite_slice, composite, sliceRegion );
+
+ // Update vectors with each image
+ data3D[ i ] = input_slice;
+ binaryTree[ i ] = labels_slice;
+ fusion[ i ] = composite_slice;
+
+ } // rof
+
+ // Random walker slice-by-slice
+ this->_GoDown( binaryTree, data3D, fusion, numSlices );
+ this->_GoUp( binaryTree, data3D, fusion, numSlices );
+
+ // Join results
+ typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
+ typename _TJoin::Pointer join = _TJoin::New( );
+ for( int i = 0; i < numSlices; ++i )
+ join->SetInput( i, binaryTree[ i ] );
+ join->Update( );
+ this->GetOutput( )->Graft( join->GetOutput( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+void
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+_Composite(
+ typename TScalarImage::Pointer& composite,
+ const TLabels* labels,
+ const TScalarImage* vesselness
+ )
+{
+ // Min-Max vesselness values
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ typename _TMinMax::Pointer minMax = _TMinMax::New( );
+ minMax->SetImage( vesselness );
+ minMax->Compute( );
+ TScalar maxVess = minMax->GetMaximum( );
+ typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
+
+ // Composite image
+ typename TScalarImage::Pointer composite = TScalarImage::New( );
+ composite->CopyInformation( vesselness );
+ composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
+ composite->Allocate( );
+ itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
+ itk::ImageRegionConstIterator< TLabels > l( labels, region );
+ itk::ImageRegionIterator< TScalarImage > c( composite, region );
+ v.GoToBegin( );
+ l.GoToBegin( );
+ c.GoToBegin( );
+ for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
+ c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+template< class _TSlicePtr, class _TInput >
+void
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+_Slice(
+ _TSlicePtr& slice,
+ const _TInput* input,
+ typename _TInput::RegionType region
+ )
+{
+ typedef typename _TSlicePTr::ObjectType _TSlice;
+ typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
+
+ typename _TExtract::Pointer extract = _TExtract::New( );
+ extract->SetInput( input );
+ extract->SetExtractionRegion( region );
+ extract->SetDirectionCollapseToIdentity( );
+ extract->Update( );
+ slice = extract->GetOutput( );
+ slice->DisconnectPipeline( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+template< class _TBinaryTree, class _TData3D, class _TFusion >
+void
+fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
+_GoDown(
+ _TBinaryTree& binaryTree,
+ const _TData3D& data3D,
+ const _TFusion& fusion,
+ int numSlices
+ )
+{
+ typedef typename _TBinaryTree::value_type _TSliceBinaryPtr;
+ typedef typename _TData3D::value_type _TSliceData3DPtr;
+ typedef typename _TFusion::value_type _TSliceFusionPtr;
+ typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
+ typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
+ typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
+
+ int z = -1;
+ while( z < numSlices - 2 )
+ {
+ z++;
+ _TSliceBinaryPtr tmp = binaryTree[ z ];
+ _TSliceBinaryPtr next = binaryTree[ z + 1];
+ _TSliceData3DPtr data = data3D[ z + 1 ];
+ _TSliceFusionPtr vess = fusion[ z + 1 ];
+
+ // Create seeds
+ _TSliceBinaryPtr seeds = _TSliceBinary::New( );
+ seeds->CopyInformation( tmp );
+ seeds->SetRequestedRegion( tmp->GetRequestedRegion( ) );
+ seeds->Allocate( );
+ seeds->FillBuffer( 0 );
+
+ /* TODO
+
+
+ SliceType::RegionType sliceRegion = tmp->GetLargestPossibleRegion();
+ SliceType::SizeType size = sliceRegion.GetSize();
+ SliceType::SpacingType spacing = tmp->GetSpacing();
+ SliceType::PointType origin = tmp->GetOrigin();
+ SliceType::IndexType start = sliceRegion.GetIndex();
+
+ SliceType::Pointer seedsObject = SliceType::New();
+ seedsObject->SetRegions(sliceRegion);
+ seedsObject->SetOrigin(origin);
+ seedsObject->SetSpacing(spacing);
+ seedsObject->Allocate();
+ seedsObject->FillBuffer(0);
+
+ SliceType::Pointer seedsBackground = SliceType::New();
+ seedsBackground->SetRegions(sliceRegion);
+ seedsBackground->SetOrigin(origin);
+ seedsBackground->SetSpacing(spacing);
+ seedsBackground->Allocate();
+ seedsBackground->FillBuffer(0);
+
+ typedef itk::ImageRegionIteratorWithIndex<SliceType> SliceIteratorType;
+ SliceIteratorType itSliceMori(tmp, tmp->GetLargestPossibleRegion());
+
+ int len = 0;
+ int numO = 0;
+ int numB = 0;
+
+ itSliceMori.GoToBegin();
+
+ while (!itSliceMori.IsAtEnd())
+ {
+ if (itSliceMori.Get() != 0)
+ len++;
+
+ typename TImage::PixelType vessVal = vess->GetPixel(itSliceMori.GetIndex());
+ typename TImage::PixelType imgVal = data->GetPixel(itSliceMori.GetIndex());
+ typename TImage::PixelType bw = tmp->GetPixel(itSliceMori.GetIndex());
+
+ if (bw != 0 && vessVal > this->m_VessTh * this->m_vessMax / 100)
+ {
+ seedsObject->SetPixel(itSliceMori.GetIndex(), 1);
+ numO++;
+ }
+
+ itSliceMori++;
+ }
+
+ if (len == 0)
+ {
+ std::cout << "Slice " << z << ": - brak pikseli zakwalifikowanych do drzewa" << std::endl;
+ continue;
+ }
+
+ SliceIteratorType itSliceVess(vess, vess->GetLargestPossibleRegion());
+ itSliceVess.GoToBegin();
+
+ while (!itSliceVess.IsAtEnd())
+ {
+ if (itSliceVess.Get() < this->m_VessTh * this->m_vessMax / 100)
+ {
+ seedsBackground->SetPixel(itSliceVess.GetIndex(), 1);
+ numB++;
+ }
+ itSliceVess++;
+
+ }
+
+ std::cout << "slice: " << z << "seeds obj: " << numO << " seeds bkg " << numB << std::endl;
+ // przygotowania do Random Walkera
+
+ typename TImage::PixelType *bufferImg = 0;
+ typename TImage::PixelType *bufferSeedsObj = 0;
+ typename TImage::PixelType *bufferSeedsBkg = 0;
+
+ int dims[2], dims2[2], dims3[2];
+ double spac[2], orig[2], spac2[2], orig2[2], spac3[2], orig3[2];
+
+ image2buffer(bufferImg, dims, spac, orig, 2, data);
+ image2buffer(bufferSeedsObj, dims2, spac2, orig2, 2, seedsObject);
+ image2buffer(bufferSeedsBkg, dims3, spac3, orig3, 2, seedsBackground);
+
+ // TO DO HERE
+ int width = dims[0];
+ int height = dims[1];
+
+ PixelType *mask = new PixelType[width*height];
+ double *probs = new double[2 * width*height];
+
+ int i, j, no_seeds = 0;
+
+ for (i = 0; i < width*height; i++)
+
+ if (bufferSeedsObj[i] == 1 || bufferSeedsBkg[i] == 1)
+ no_seeds++;
+
+ int *seed_indexes = new int[no_seeds];
+ int *seed_labels = new int[no_seeds];
+
+ for (j = 0, i = 0; i < width*height; i++)
+ {
+ if (bufferSeedsObj[i] == 1)
+ {
+ seed_indexes[j] = i + 1;//=i+1 (i)
+ seed_labels[j] = this->m_replaceValue;
+ j++;
+ }
+ else if (bufferSeedsBkg[i] == 1)
+ {
+ seed_indexes[j] = i + 1;//=i+1 (i)
+ seed_labels[j] = 0;
+ j++;
+ }
+ }
+
+ call_walker(mask, probs, bufferImg, height, width, seed_indexes, no_seeds, seed_labels, 2, this->m_beta);
+
+ delete[] bufferSeedsObj;
+ delete[] bufferSeedsBkg;
+ delete[] seed_indexes;
+ delete[] seed_labels;
+
+ SliceType::Pointer img2 = buffer2image(mask, dims, spac, orig, 2);
+
+ SliceType::PixelType *probs_img = new PixelType[width*height];
+ for (i = 0; i < width*height; i++)
+ probs_img[i] = floor(255 * probs[i]);
+
+ SliceType::Pointer img3 = buffer2image(probs_img, dims, spac, orig, 2);
+
+ typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
+ MaxFilterType::Pointer maxFilter = MaxFilterType::New();
+ maxFilter->SetInput(0, img2);
+ maxFilter->SetInput(1, next);
+ maxFilter->Update();
+
+ SliceType::Pointer sth = maxFilter->GetOutput();
+ ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());
+
+ delete[] probs;
+ */
+ }
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalarImage >
+template< class _TBinaryTree, class _TData3D, class _TFusion >
+void
+_GoUp(
+ _TBinaryTree& binaryTree,
+ const _TData3D& data3D,
+ const _TFusion& fusion,
+ int numSlices
+ )
+{
+}
+
+#endif // __fpa__Common__SliceBySliceRandomWalker__hxx__
+// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff