[Insight-users] Fwd: Re: Filling big 3D holes

[Gib Bogle]

(Sorry about forgetting to "reply-all")

-------- Original Message --------
Subject: Re: [Insight-users] Filling big 3D holes
Date: Thu, 03 Feb 2011 17:06:42 +1300
From: Gib Bogle <g.bogle at auckland.ac.nz>
To: brian avants <stnava at gmail.com>

Hi Brian,

Before I delve into this, do you think it overcomes the problem of wall
continuity that I refer to in my reply to Richard?

Thanks
Gib

On 3/02/2011 5:00 p.m., brian avants wrote:
> hi gib
>
> following on richard's response
>
>    // hole-filling algorithm :
>     1. get distance map of background from object
>     2. threshold above zero
>     3. label connected components
>     4. label surface of each connected component
>     5. if everywhere on the surface is next to object then it's a hole
>
> an itk code sketch that approximates this is below ---
>
> brian
>
>
> template<unsigned int ImageDimension>
> int FillHoles(int argc, char *argv[])
> {
>
>    typedef itk::DanielssonDistanceMapImageFilter<ImageType, ImageType>
>   FilterType;
>    typename  FilterType::Pointer filter = FilterType::New();
>    filter->InputIsBinaryOff();
>    filter->SetUseImageSpacing(false);
>    filter->SetInput(image1);
>    filter->Update();
>
>    typedef itk::Image<int, ImageDimension>  labelimagetype;
>    typename ImageType::Pointer dist = filter->GetOutput();
>    typename ImageType::Pointer regions =
> BinaryThreshold<ImageType>(0.001,1.e9,1,dist);
>
>    typedef itk::ConnectedComponentImageFilter<  labelimagetype,
> labelimagetype>  ccFilterType;
>    typedef itk::RelabelComponentImageFilter<  labelimagetype, ImageType
>> RelabelType;
>    typename RelabelType::Pointer relabel = RelabelType::New();
>    typename ccFilterType::Pointer ccfilter = ccFilterType::New();
>
>    typename CastFilterType::Pointer castRegions = CastFilterType::New();
>    castRegions->SetInput( regions );
>
>    ccfilter->SetInput(castRegions->GetOutput());
>    ccfilter->SetFullyConnected( 0 );
>    ccfilter->Update();
>    relabel->SetInput( ccfilter->GetOutput() );
>
>    relabel->SetMinimumObjectSize( 0 );
>    //    relabel->SetUseHistograms(true);
>    try
>      {
>      relabel->Update();
>      }
>    catch( itk::ExceptionObject&  excep )
>      {
>      std::cerr<<  "Relabel: exception caught !"<<  std::endl;
>      std::cerr<<  excep<<  std::endl;
>      }
>
>
>      // WriteImage<ImageType>(relabel->GetOutput(),"test.nii");
>
>    if (holeparam == 2)
>      {
>        std::cout<<" Filling all holes"<<   std::endl;
>        typedef itk::ImageRegionIteratorWithIndex<ImageType>  Iterator;
>        Iterator vfIter( relabel->GetOutput(),
> relabel->GetOutput()->GetLargestPossibleRegion() );
>
>        for(  vfIter.GoToBegin(); !vfIter.IsAtEnd(); ++vfIter )
> 	{
> 	  if (vfIter.Get()>  1 ) imageout->SetPixel(vfIter.GetIndex(),1);
> 	}
>
> //      WriteImage<ImageType>(imageout,outname.c_str());
>
>        return 0;
>      }
>
>
>    typedef itk::NeighborhoodIterator<ImageType>   iteratorType;
>    typename iteratorType::RadiusType rad;
>    for (unsigned int j=0; j<ImageDimension; j++) rad[j]=1;
>    iteratorType GHood(rad,
> relabel->GetOutput(),relabel->GetOutput()->GetLargestPossibleRegion());
>
>    float maximum=relabel->GetNumberOfObjects();
>    //now we have the exact number of objects labeled independently
>    for (unsigned int lab=2; lab<=maximum; lab++)
>      {
>        float erat=2;
>        if (holeparam<= 1 )
> 	{
>        GHood.GoToBegin();
>
>        unsigned long objectedge=0;
>        unsigned long backgroundedge=0;
>        unsigned long totaledge=0;
>        unsigned long volume=0;
>
>        while (!GHood.IsAtEnd())
> 	{
> 	  typename ImageType::PixelType p = GHood.GetCenterPixel();
> 	  typename ImageType::IndexType ind = GHood.GetIndex();
> 	  typename ImageType::IndexType ind2;
> 	  if ( p == lab )
> 	    {
> 	      volume++;
> 	      for (unsigned int i = 0; i<  GHood.Size(); i++)
> 		{
> 		  ind2=GHood.GetIndex(i);
> 		  float dist=0.0;
> 		  for (unsigned int j=0; j<ImageDimension; j++)
> 		    dist+=(float)(ind[j]-ind2[j])*(float)(ind[j]-ind2[j]);
> 		  dist=sqrt(dist);
> 		  float val2=image1->GetPixel(ind2);
> 		  if (val2>= 0.5&&  GHood.GetPixel(i) != lab )
> 		    {
> 		      objectedge++;
> 		      totaledge++;
> 		    }
> 		  else if (val2<  0.5&&  GHood.GetPixel(i) != lab )
> 		    {
> 		      backgroundedge++;
> 		      totaledge++;
> 		    }
> 		}
> 	    }
> 	  ++GHood;
> 	}
>        float vrat=(float)totaledge/(float)volume;
>        erat=(float)objectedge/(float)totaledge;
>        std::cout<<" Lab "<<  lab<<  " volume "<<  volume<<  " v-rat"
> <<  vrat<<  " edge "<<  erat<<  std::endl;
> 	}
>
>        if (erat>  holeparam)// fill the hole
> 	{
> 	  std::cout<<" Filling "<<  lab<<  " of "<<  maximum<<   std::endl;
> 	  typedef itk::ImageRegionIteratorWithIndex<ImageType>  Iterator;
> 	  Iterator vfIter( relabel->GetOutput(),
> relabel->GetOutput()->GetLargestPossibleRegion() );
>
> 	  for(  vfIter.GoToBegin(); !vfIter.IsAtEnd(); ++vfIter )
> 	    {
> 	      if (vfIter.Get() == lab ) imageout->SetPixel(vfIter.GetIndex(),1);
> 	    }
>
> 	}
>
>      }
>
> // write the image
> //  WriteImage<ImageType>(imageout,outname.c_str());
>
>
>    return 0;
>
> }
>
>
>
>
>
>
> On Wed, Feb 2, 2011 at 10:28 PM, Gib Bogle<g.bogle at auckland.ac.nz>  wrote:
>> I'm working with a volume image that was generated by labelling the laminae
>> of blood vessels.  My aim is to segment out the vasculature.  There are many
>> difficulties, and the particular issue I'm addressing at the moment is
>> filling in the vessels.  The intensity of the labelling of the walls is
>> variable, with patches that are indistinguishable from background.  The
>> vessel diameters vary widely, from about 4 to about 60 voxels.  After some
>> preprocessing I have a binary image, on which the ITK hole-filling function
>> works well with the small-diameter vessels, but the big vessels present a
>> problem, even when the walls are "watertight" (i.e. without holes).
>>
>> My best idea so far is to send probes out in all 26 directions (all
>> neighbours of a voxel) and count the number of probes that hit a wall within
>> a specified radius.  It's tricky to specify both the radius and the critical
>> number of hits, without getting too many false positives (voxels outside the
>> vessels showing up as inside).  (The filter is of course applied
>> iteratively.)
>>
>> I'm wondering if anyone else here has addressed a similar problem.
>>
>> Thanks
>> Gib
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.html
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>