[Insight-users] Filling big 3D holes

[brian avants]

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
>