#define __fpa__Common__SliceBySliceRandomWalker__hxx__
#include <vector>
+#include <itkBinaryThresholdImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkJoinSeriesImageFilter.h>
+#include <itkMaximumImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
+#include <fpa/Common/RandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
SliceBySliceRandomWalker( )
- : Superclass( )
+ : Superclass( ),
+ m_Epsilon( TScalar( 1e-5 ) ),
+ m_Beta( TScalar( 1 ) ),
+ m_VesselnessThreshold( TScalar( 5 ) )
{
ivqITKInputConfigureMacro( InputLabels, TLabels );
ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
const TScalarImage* vesselness = this->GetInputVesselness( );
this->AllocateOutputs( );
+ // Get vesselness range
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ typename _TMinMax::Pointer vMinMax = _TMinMax::New( );
+ vMinMax->SetImage( vesselness );
+ vMinMax->Compute( );
+ TScalar vMax = vMinMax->GetMaximum( );
+
// Some values
typename TImage::RegionType region = input->GetRequestedRegion( );
- typename TImage::SizeType regionIndex = region.GetIndex( );
+ typename TImage::IndexType regionIndex = region.GetIndex( );
typename TImage::SizeType sliceSize = region.GetSize( );
int numSlices = sliceSize[ 2 ];
sliceSize[ 2 ] = 0;
+ this->m_VesselnessValue = TScalar( 0.01 );
+ this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax;
// Composite image
- typename TScalarImage::Pointer composite = TScalarImage::New( );
- composite->CopyInformation( vesselness );
- composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
- composite->Allocate( );
- this->_Composite( composite, labels, vesselness );
+ typename TScalarImage::Pointer composite;
+ this->_Composite( composite, labels, vesselness, vMax );
// Extract slices
std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
typename _TSliceScalarImage::Pointer composite_slice;
this->_Slice( input_slice, input, sliceRegion );
this->_Slice( labels_slice, labels, sliceRegion );
- this->_Slice( composite_slice, composite, sliceRegion );
+ this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion );
// Update vectors with each image
data3D[ i ] = input_slice;
} // rof
// Random walker slice-by-slice
- this->_GoDown( binaryTree, data3D, fusion, numSlices );
- this->_GoUp( binaryTree, data3D, fusion, numSlices );
+ this->_RandomWalker( binaryTree, data3D, fusion, true );
+ this->_RandomWalker( binaryTree, data3D, fusion, false );
// Join results
- typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
+ typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin;
typename _TJoin::Pointer join = _TJoin::New( );
for( int i = 0; i < numSlices; ++i )
join->SetInput( i, binaryTree[ i ] );
+ join->SetOrigin( input->GetOrigin( )[ 2 ] );
+ join->SetSpacing( input->GetSpacing( )[ 2 ] );
join->Update( );
this->GetOutput( )->Graft( join->GetOutput( ) );
}
_Composite(
typename TScalarImage::Pointer& composite,
const TLabels* labels,
- const TScalarImage* vesselness
+ const TScalarImage* vesselness,
+ const TScalar& maxVess
)
{
- // Min-Max vesselness values
- typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
- typename _TMinMax::Pointer minMax = _TMinMax::New( );
- minMax->SetImage( vesselness );
- minMax->Compute( );
- TScalar maxVess = minMax->GetMaximum( );
+ // Some values
typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
// Composite image
- typename TScalarImage::Pointer composite = TScalarImage::New( );
+ composite = TScalarImage::New( );
composite->CopyInformation( vesselness );
composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
composite->Allocate( );
typename _TInput::RegionType region
)
{
- typedef typename _TSlicePTr::ObjectType _TSlice;
+ typedef typename _TSlicePtr::ObjectType _TSlice;
typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
typename _TExtract::Pointer extract = _TExtract::New( );
template< class _TBinaryTree, class _TData3D, class _TFusion >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
-_GoDown(
+_RandomWalker(
_TBinaryTree& binaryTree,
const _TData3D& data3D,
const _TFusion& fusion,
- int numSlices
+ bool down
)
{
typedef typename _TBinaryTree::value_type _TSliceBinaryPtr;
typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
+ typedef typename _TSliceFusion::PixelType _TSliceScalar;
- int z = -1;
- while( z < numSlices - 2 )
+ int numSlices = binaryTree.size( );
+ int z = ( down )? -1: numSlices;
+ int o = ( down )? 1: -1;
+ while( ( down && z < numSlices - 2 ) || ( !down && z > 1 ) )
{
- z++;
+ z += o;
_TSliceBinaryPtr tmp = binaryTree[ z ];
- _TSliceBinaryPtr next = binaryTree[ z + 1];
- _TSliceData3DPtr data = data3D[ z + 1 ];
- _TSliceFusionPtr vess = fusion[ z + 1 ];
+ _TSliceBinaryPtr next = binaryTree[ z + o ];
+ _TSliceData3DPtr data = data3D[ z + o ];
+ _TSliceFusionPtr vess = fusion[ z + o ];
+ typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
- // Create seeds
+ // Create seeds and background
_TSliceBinaryPtr seeds = _TSliceBinary::New( );
seeds->CopyInformation( tmp );
- seeds->SetRequestedRegion( tmp->GetRequestedRegion( ) );
+ seeds->SetBufferedRegion( region );
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
- )
-{
+ itk::ImageRegionConstIterator< _TSliceBinary > tIt( tmp, region );
+ itk::ImageRegionConstIterator< _TSliceFusion > vIt( vess, region );
+ itk::ImageRegionIterator< _TSliceBinary > sIt( seeds, region );
+ tIt.GoToBegin( );
+ vIt.GoToBegin( );
+ sIt.GoToBegin( );
+ int len = 0;
+ int numO = 0;
+ int numB = 0;
+ for( ; !tIt.IsAtEnd( ); ++tIt, ++vIt, ++sIt )
+ {
+ if( tIt.Get( ) != 0 )
+ {
+ len++;
+ if( vIt.Get( ) > this->m_VesselnessValue )
+ {
+ sIt.Set( 1 );
+ numO++;
+
+ } // fi
+
+ } // fi
+
+ if( vIt.Get( ) < this->m_VesselnessValue )
+ {
+ sIt.Set( 2 );
+ numB++;
+
+ } // fi
+
+ } // rof
+
+ if( len == 0 )
+ continue;
+
+ // Random walker function
+ typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
+ typename _TFunction::Pointer edge = _TFunction::New( );
+ edge->SetBeta( this->m_Beta );
+ edge->SetEpsilon( this->m_Epsilon );
+ edge->TreatAsWeightOff( );
+
+ // Real random walker
+ typedef fpa::Common::RandomWalker< _TSliceData3D, _TSliceBinary, _TSliceScalar > _TRandomWalker;
+ typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
+ rw->SetInput( data );
+ rw->SetInputLabels( seeds );
+ rw->SetEdgeFunction( edge );
+
+ // Extract label
+ typedef itk::BinaryThresholdImageFilter< _TSliceBinary, _TSliceBinary > _TExtract;
+ typename _TExtract::Pointer extract = _TExtract::New( );
+ extract->SetInput( rw->GetOutput( ) );
+ extract->SetInsideValue( 1 );
+ extract->SetOutsideValue( 0 );
+ extract->SetLowerThreshold( 1 );
+ extract->SetUpperThreshold( 1 );
+
+ // Keep maximum
+ typedef itk::MaximumImageFilter< _TSliceBinary > _TMax;
+ typename _TMax::Pointer max = _TMax::New();
+ max->SetInput( 0, extract->GetOutput( ) );
+ max->SetInput( 1, next );
+ max->Update( );
+ binaryTree[ z + o ] = max->GetOutput( );
+ binaryTree[ z + o ]->DisconnectPipeline( );
+
+ } // elihw
}
#endif // __fpa__Common__SliceBySliceRandomWalker__hxx__
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff