X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FCommon%2FSliceBySliceRandomWalker.hxx;h=3b86d023b9b56f65de65e71121648c6fafce3025;hb=463b1ec45c70ca63aa9f19f8cb58ae5d5134e56b;hp=54b1fbcfcebdc2735270a05d7754dc3b97fd0f98;hpb=0251081f09de5bd702c01565c9401c1eb1983ae5;p=FrontAlgorithms.git diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.hxx b/lib/fpa/Common/SliceBySliceRandomWalker.hxx index 54b1fbc..3b86d02 100644 --- a/lib/fpa/Common/SliceBySliceRandomWalker.hxx +++ b/lib/fpa/Common/SliceBySliceRandomWalker.hxx @@ -6,17 +6,24 @@ #define __fpa__Common__SliceBySliceRandomWalker__hxx__ #include +#include #include #include #include #include +#include #include +#include +#include // ------------------------------------------------------------------------- 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 ); @@ -45,19 +52,25 @@ GenerateData( ) 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 ); @@ -76,7 +89,7 @@ GenerateData( ) 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; @@ -86,14 +99,16 @@ GenerateData( ) } // 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( ) ); } @@ -105,19 +120,15 @@ fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >:: _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( ); @@ -142,7 +153,7 @@ _Slice( 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( ); @@ -159,11 +170,11 @@ template< class _TImage, class _TLabels, class _TScalarImage > 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; @@ -172,180 +183,95 @@ _GoDown( 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 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 - ) -{ + 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__