[Insight-users] Skeleton Filter?

[andreas.krebs at rad.uni-kiel.de]

Hi Zihua Su,


Zihua Su schrieb:
> > Hi Andreas,
> > Thank you very much for your detailed and kind help .
> > Sorry , i didn't make myself clear. I am processing a data 512*512*352
> > sclice by slice. I have tried BinaryThinningImageFilter which is
> > straight forward to find the skeleton in an images. The problem is
> > that it just takes too long.
> > For one slice, it takes approximately 10mins on my computer.
This really sounds too slow. But most skeletonization algorithm work
with a queue which becomes to big when the image size grows. This seems
to be the problem when I use Lamy's Digital Topology algorithm for
512x512x400 on a 8 GByte RAM machine. So, I splitted the whole image
into four boxes with 100 sclices. The skeletonization of 512x512x100
costs about 4500 seconds of CPU time and needed 6.5 GByte. Of course,
this is not the same as a skeletonization of the whole image. And
anyway, it is too expensive.

But i
> > have also tried chamfer distance in ITK which gave me a beautiful
> > distance map for the whole 3D data resonable quickly. So, i think my
> > goal now will be take out the skeleton from the distance map.

> > You euqtion about "for all q \in N_26(p): d(p) \neq d(q) - Dist(p,q)"
> > is interesting for me, looks like what i want. Could you pls give me a
> > reference for this. I guess it is calculating some vablue from the
> > relation of neighbour points.

The reference is the mentioned paper of Chris Pudney, p. 408.
Furthermore, the following paper may help for a justification of the
maximal ball approach:

T. Hildebrand, P. Rüegsegger: A new method for the model-independen
assessment of thickness in three-dimensional images.

I attach a code sequence. I think the fuzzy distance transformation
image can be easily replaced by your distance transform.

And ( mbsImIt.Get() + mbsIm->GetPixel( nb ) ) by 2.

hdeltad should be replaced by your chamfer weights, probably 3,4,5 if
you are using isotropic voxels.

The idea of the maximal balls is simple. Just look at all your
neighbors. If they are smaller your are a center of a max. ball.

Bests Andreas

> >
> > Regards
> >
> > Zihua Su
> >
> > On 4/6/07, Andreas Krebs <krebs.andreas at gmx.de> wrote:
>> >> Hi Zihua Su
>> >>
>> >> I already tried the filter proposed in the document
>> >> Digital Topology from Julien Lamy. It works fine, but the
>> >> terminality criterion needed to be improved for my application. I
did it
>> >> and it works, but I programmed it in a very ugly way. So all in all
>> >> I have an ugly implementation of
>> >>
>> >> Chris Pudney: "Distance-Ordered Homotopic Thinning:
>> >> A Skeletonization Algorithm for 3D Digital Images"
>> >> in: Computer Vision and Image Understanding, Vol 72, No. 3,
>> >> December, pp. 404--413, 1998
>> >>
>> >> which works fine for images of 100 x 100 x 100. It becomes very slow
for
>> >> bigger images.
>> >>
>> >> If you already calculate the ChamferDistanceMap, then it is  worth to
>> >> ask, if you really need the skeleton or if the skeleton like ridges
>> >> which can be calculated (and implemented) much faster are sufficient.
>> >>
>> >> For my application it was sufficient to look for the centers of maximal
>> >> disks (spheres). These can be identified from the chamfer distance
>> >> transform d(.) . The centers are the points p that satisfy
>> >>
>> >> for all q \in N_26(p): d(p) \neq d(q) - Dist(p,q)
>> >>
>> >> where Dist(p,q) is the distance between p and q (cf. Pudney, p. 408,
eq.
>> >> (3) ). I use this approach also to find ridges in an grayscale image
>> >> based on fuzzy distance transform fdt(.). The output looks almost
like a
>> >> skeleton.
>> >>
>> >> If somebody would like to put these routines into the itk template
>> >> filter design, I would be glad. Perhaps my code can be helpful to
>> >> fulfill this task.
>> >>
>> >> Furthermore, there is given a more complicated method to identify the
>> >> ridge using the Laplacian somewhere in the InsightApplicatiosn - I can
>> >> not find the link now. Anyway, I prefer my fuzzy distance
transformation
>> >> ridges, after I played around with the Laplacian method.
>> >>
>> >> Regards
>> >>        Andreas
>> >>
>> >> Zihua Su schrieb:
>>> >>> Hi all,
>>> >>>
>>> >>> I have implemented the fastchamferdistancemap sucessfully. but i am
>>> >>> wondering if there is a filter to take out the skeleton out in ITK!
>>> >>> Any clue will be welcome and appreciated.
>>> >>>
>>> >>> Cheers
>>> >>>



float hdeltad[6]; // for hdelta in fdt

inline float hdelta (int i, int j, int k) {
  return hdeltad[(abs(i) +abs(j)) +3*abs(k)];
}



void Inithdelta ( const ImageType::SpacingType sp) {
  const float hx=sp[0], hy=sp[1], hz=sp[2];

  hdeltad[0]= 0.;
  hdeltad[1]= hx/2.;
  hdeltad[2]= sqrt(2*hx*hx)/2.;
  hdeltad[3]= hz/2.;
  hdeltad[4]= sqrt(hx*hx + hz*hz)/2.;
  hdeltad[5]= sqrt(2*hx*hx + hz*hz)/2.;

  return;
}

BinaryImageType::Pointer centerofball( ImageType::Pointer mbsIm,
				       ImageType::Pointer iIm, float thresh) {
  //mbsIm = original grayscale image
  //iIm = fuzzy distance transform of mbsIm
  //thresh = threshold for minimum ball size relative to voxelsize in x
direction
  BinaryImageType::Pointer oIm = BinaryImageType::New();
  BinaryImageType::RegionType oRegion;
  std::cout<<"Compute centers of balls"<<std::endl;
  oRegion = iIm->GetLargestPossibleRegion();
  oIm->SetRegions( oRegion );
  oIm->SetSpacing( iIm->GetSpacing() );
  oIm->SetOrigin( iIm->GetOrigin() );
  oIm->Allocate();

  const ImageType::RegionType::SizeType& size = oRegion.GetSize();
  itk::ImageRegionConstIteratorWithIndex< ImageType> mbsImIt( mbsIm,
oRegion );
  itk::ImageRegionConstIteratorWithIndex< ImageType> iImIt( iIm, oRegion );
  itk::ImageRegionIterator< BinaryImageType> oImIt( oIm, oRegion );

  ImageType::IndexType index, nb;

  for ( iImIt.GoToBegin(), mbsImIt.GoToBegin(), oImIt.GoToBegin();
	!iImIt.IsAtEnd();
	++iImIt, ++mbsImIt, ++oImIt) {
    index = iImIt.GetIndex();
    float fdtvalue = iImIt.Get();
    float mbvalue = mbsImIt.Get();

    // cf. Eq.(3) from p. 408, Chris Pudney, Homotopic Thinning
    bool center = fdtvalue > thresh *hdeltad[1]; // ball is sufficiently big
    // loop over all neighbors until center is denied
    for (nb[2]=index[2]-1; center && nb[2]<=index[2]+1; nb[2]++) {
      if ( nb[2] < 0 | nb[2] >= size[2] ) continue;

      for (nb[1]=index[1]-1; center && nb[1]<=index[1]+1; nb[1]++) {
	if ( nb[1] < 0 | nb[1] >= size[1] ) continue;

	for (nb[0]=index[0]-1; center && nb[0]<=index[0]+1; nb[0]++) {
	  if ( nb[0] < 0 | nb[0] >= size[0] ) continue;
	  if ( nb[2]==index[2] && nb[1]==index[1] && nb[0]==index[0] )
	    continue; //=> next neighbour

	  if( iIm->GetPixel( nb ) > fdtvalue
	      + 0.6 *( mbsImIt.Get() + mbsIm->GetPixel( nb ) )
	      * hdelta( nb[0]-index[0],
			nb[1]-index[1],
			nb[2]-index[2]) )
	    center = false; //not a center point
	}
      }
    }
    if ( center )
      oImIt.Set( 255 ); // center of a ball
    else
      oImIt.Set( 0 );
  }
  return oIm;
}