From 03a2bd7a6ab76afddd40da4c85bfeb7869b2c410 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Thu, 12 Oct 2017 21:19:48 -0500 Subject: [PATCH] ... --- lib/fpa/Common/SliceBySliceRandomWalker.hxx | 194 ++++++-------------- 1 file changed, 59 insertions(+), 135 deletions(-) diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.hxx b/lib/fpa/Common/SliceBySliceRandomWalker.hxx index 54b1fbc..d1fbd9c 100644 --- a/lib/fpa/Common/SliceBySliceRandomWalker.hxx +++ b/lib/fpa/Common/SliceBySliceRandomWalker.hxx @@ -11,6 +11,8 @@ #include #include #include +#include +#include // ------------------------------------------------------------------------- template< class _TImage, class _TLabels, class _TScalarImage > @@ -51,6 +53,8 @@ GenerateData( ) typename TImage::SizeType sliceSize = region.GetSize( ); int numSlices = sliceSize[ 2 ]; sliceSize[ 2 ] = 0; + this->m_VesselnessValue = + this->m_VesselnessThreshold * this->m_VesselnessMax / TScalar( 100 ); // Composite image typename TScalarImage::Pointer composite = TScalarImage::New( ); @@ -172,155 +176,76 @@ _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 ) { z++; _TSliceBinaryPtr tmp = binaryTree[ z ]; - _TSliceBinaryPtr next = binaryTree[ z + 1]; + _TSliceBinaryPtr next = binaryTree[ z + 1 ]; _TSliceData3DPtr data = data3D[ z + 1 ]; _TSliceFusionPtr vess = fusion[ z + 1 ]; + typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( ); - // Create seeds + // Create seeds and background _TSliceBinaryPtr seeds = _TSliceBinary::New( ); seeds->CopyInformation( tmp ); - seeds->SetRequestedRegion( tmp->GetRequestedRegion( ) ); + seeds->SetRequestedRegion( region ); seeds->Allocate( ); seeds->FillBuffer( 0 ); + 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_Eps ); + 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 ); + + // Keep maximum /* 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); @@ -329,7 +254,6 @@ _GoDown( SliceType::Pointer sth = maxFilter->GetOutput(); ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion()); - delete[] probs; */ } -- 2.45.1