From 0251081f09de5bd702c01565c9401c1eb1983ae5 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Thu, 12 Oct 2017 16:34:26 -0500 Subject: [PATCH] ... --- appli/CTBronchi/FastRandomWalker.cxx | 50 ++- lib/fpa/Common/SliceBySliceRandomWalker.h | 93 ++++++ lib/fpa/Common/SliceBySliceRandomWalker.hxx | 352 ++++++++++++++++++++ 3 files changed, 467 insertions(+), 28 deletions(-) create mode 100644 lib/fpa/Common/SliceBySliceRandomWalker.h create mode 100644 lib/fpa/Common/SliceBySliceRandomWalker.hxx diff --git a/appli/CTBronchi/FastRandomWalker.cxx b/appli/CTBronchi/FastRandomWalker.cxx index b2c901b..7342562 100644 --- a/appli/CTBronchi/FastRandomWalker.cxx +++ b/appli/CTBronchi/FastRandomWalker.cxx @@ -6,7 +6,8 @@ #include #include #include -#include "MoriLabelling.h" +#include +#include #include "Functions.h" // ------------------------------------------------------------------------- @@ -16,29 +17,25 @@ 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; // ------------------------------------------------------------------------- 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 ); @@ -65,26 +62,23 @@ int main( int argc, char* 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 ) { diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.h b/lib/fpa/Common/SliceBySliceRandomWalker.h new file mode 100644 index 0000000..48f2bed --- /dev/null +++ b/lib/fpa/Common/SliceBySliceRandomWalker.h @@ -0,0 +1,93 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __fpa__Common__SliceBySliceRandomWalker__h__ +#define __fpa__Common__SliceBySliceRandomWalker__h__ + +#include +#include + +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 +#endif // ITK_MANUAL_INSTANTIATION +#endif // __fpa__Common__SliceBySliceRandomWalker__h__ +// eof - $RCSfile$ diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.hxx b/lib/fpa/Common/SliceBySliceRandomWalker.hxx new file mode 100644 index 0000000..54b1fbc --- /dev/null +++ b/lib/fpa/Common/SliceBySliceRandomWalker.hxx @@ -0,0 +1,352 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__ +#define __fpa__Common__SliceBySliceRandomWalker__hxx__ + +#include +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +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 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 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$ -- 2.45.1