[Insight-users] Skeleton Filter?
Andreas Krebs, UKSH
andreas.krebs at rad.uni-kiel.de
Tue Apr 10 16:01:17 EDT 2007
Dear Ivan Macia,
in the following I try my best to share it - the code is not ready to be
packed, so I just copy the necessary routines of source code, hope this
helps.
Concerning the ridge extraction: Firstly, you need a chamfer distance
transformation (e.g. that given in Lamy's Digital Topology paket) to get
a distance map of the object, you would like to extract the ridge from.
If you are working with fuzzy distance objects, I can send you also a
fuzzy distance algorithm of
> @article{SaChMa97,
> author = "Saha, Punam K. and Chaudhuri, B. B. and Dutta
> Majumder, D.",
> title = "A new shape preserving parallel thinning algorithm
> for 3D digital images",
> language = "English",
> journal = "Pattern Recognition",
> volume = "30",
> number = "12",
> pages = "1939-1955",
> year = {1997},
> keywords = "{3D digital topology; simple point; outer-layer;
> shape-point; skeleton; arc-skeleton; 3D parallel
> thinning; sub-fields; shape distance}",
I have not controlled it, but the chamfer distance should be quicker
than my implementation.
Secondly, you put the binary, the distance map, and a threshold into the
subroutine
BinaryImageType::Pointer centerofball( ImageType::Pointer mbsIm,
ImageType::Pointer iIm, float thresh)
I already mailed the code of the subroutine on 04/06, 21h38, and a small
description of the most important stuff to be changed. ImageType was a
FLOAT image in my case, BinaryImageType of type CHAR.
Concerning the skeleton is becomes more complicated since I have not
repaired the last bugs and I have no time to do it the next two weeks.
But maybe you can catch the idea from the following code and Chris
Pudney's publication. I use the opportunity in Lamy's code to change the
terminality criterion. I attach the Code for the two new terminality
criteria NonWitnessTerminality and SurfaceTerminality.
Perhaps the bug is in the implementation of the
NonWitnessTerminalityCriterion
// cf. Eq.(3) from p. 408, C. Pudney,
if( fabs( mod_nonwitnessorder->GetPixel(index+offset)
- hdelta(offset[0], offset[1], offset[2])
- fdtvalue ) < delta )
Also the Modul variables of the NonWitnessTerminalityCrit. have
to bypass Lamy's routines which seems dirty style to me, but it was my
quick workaround:
//Modul variables and functions
float hdeltad[6]; // for hdelta checking nonwitness points
float mod_witnessthresh;
The SurfaceTerminality criterion should work fine.
Here the calling code sequence:
typedef itk::SkeletonizeImageFilter
<BinaryImageType, itk::Connectivity<Dimension, 0> > Skeletonizer;
Skeletonizer::Pointer skeletonizer = Skeletonizer::New();
Skeletonizer::Pointer skeletonizer2 = Skeletonizer::New();
itk::ImageFunction<BinaryImageType, bool >::Pointer
nonwitnessterminality;
nonwitnessterminality = itk::NonwitnessTerminalityImageFunction
<BinaryImageType, itk::Connectivity<Dimension, 0> >::New();
itk::ImageFunction<BinaryImageType, bool >::Pointer surfaceterminality;
surfaceterminality = itk::SurfaceTerminalityImageFunction
<BinaryImageType, itk::Connectivity<Dimension, 0> >::New();
typedef itk::ChamferDistanceTransformImageFilter
< BinaryImageType, Skeletonizer::OrderingImageType>
DistanceMapFilterType;
DistanceMapFilterType::Pointer distanceMapFilter
= DistanceMapFilterType::New();
unsigned int weights[] = { 3, 4, 5 };
distanceMapFilter->SetDistanceFromObject(false);
distanceMapFilter->SetWeights(weights, weights+3);
distanceMapFilter->SetInput(binthresh->GetOutput());
distanceMapFilter->Update();
skeletonizer->SetOrderingImage(distanceMapFilter->GetOutput());
skeletonizer->SetOrderingImage(rescale->GetOutput());
/* Distance-Ordered Homotopic Thinning
cf. 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
*/
// Stage 1, Rule B with threshhold tau = skeletonstage1_arg *
hdeltad[1]
bool didstage1=false;
if ( GetInputArgs()->skeletonstage1_arg >= 0.) {
std::cout
<<"Stage1 of Distance Ordered Homotopic Thinning with threshold: "
<<GetInputArgs()->skeletonstage1_arg
<<std::endl;
skeletonizer->SetTerminalityCriterion( nonwitnessterminality );
itk::SetNonwitnessOrdering( iIm, GetInputArgs()
->skeletonstage1_arg*hdeltad[1] );
nonwitnessterminality->SetInputImage( binthresh->GetOutput() );
skeletonizer->Update();
didstage1=true;
oIm =skeletonizer->GetOutput();
}
// Stage 2, Rule A for medial surface
if ( stage2_arg == std::string("medsurface") ) {
std::cout <<"Stage2 of Distance Ordered Homotopic Thinning: "
<<stage2_arg <<std::endl;
skeletonizer2->SetTerminalityCriterion( surfaceterminality );
skeletonizer2->SetOrderingImage(rescale->GetOutput());
if ( didstage1 ) {
skeletonizer2->SetInput(skeletonizer->GetOutput() );
surfaceterminality->SetInputImage( skeletonizer->GetOutput() );
}
else {
skeletonizer2->SetInput( binthresh->GetOutput() );
surfaceterminality->SetInputImage( binthresh->GetOutput() );
}
skeletonizer2->Update();
oIm = skeletonizer2->GetOutput();
}
// Hope this helps, regards
Andreas
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NonwitnessTerminality.h
Type: text/x-chdr
Size: 2064 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20070410/342124d1/NonwitnessTerminality.h
-------------- next part --------------
/* Distance-Ordered Homotopic Thinning
cf. 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
The nonwitness terminality criterion implements Eq.(3) from p. 408
*/
#ifndef NonwitnessTerminality_txx
#define NonwitnessTerminality_txx
#include "NonwitnessTerminality.h"
namespace itk
{
//Modul variables and functions
float hdeltad[6]; // for hdelta checking nonwitness points
float mod_witnessthresh;
ImageType::Pointer mod_nonwitnessorder;
void SetNonwitnessOrdering(ImageType::Pointer const iIm, float const thresh) {
mod_nonwitnessorder = iIm;
mod_witnessthresh = thresh;
ImageType::SpacingType sp = iIm->GetSpacing();
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.;
}
inline float hdelta (int i, int j, int k) {
return hdeltad[(abs(i) +abs(j)) +3*abs(k)];
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
NonwitnessTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::NonwitnessTerminalityImageFunction()
{
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
NonwitnessTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::Evaluate(PointType const & point) const
{
typename TImage::IndexType index;
ConvertPointToNearestIndex(point, index);
return EvaluateAtIndex(index);
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
NonwitnessTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::EvaluateAtIndex(IndexType const & index) const {
TForegroundConnectivity const & fgc = TForegroundConnectivity::GetInstance();
// int nbNeighbors = 0;
float delta = hdeltad[1] /3.;
float fdtvalue = mod_nonwitnessorder->GetPixel(index);
if ( fdtvalue <= mod_witnessthresh )
return false; // => deletable
for(int i=0; i<fgc.GetNumberOfNeighbors(); ++i) {
itk::Offset<TForegroundConnectivity::Dimension> offset;
for(unsigned int j=0; j<TForegroundConnectivity::Dimension; ++j)
{
offset[j] = fgc.GetNeighborsPoints()[i][j];
}
// cf. Eq.(3) from p. 408
if( fabs( mod_nonwitnessorder->GetPixel(index+offset)
- hdelta(offset[0], offset[1], offset[2])
- fdtvalue ) < delta )
return false; // => deletable, since not center of a ball
}
// => non deletable, since center of a ball (nonwitness
//and fdtvalue > mod_witnessthresh
return true;
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
NonwitnessTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::EvaluateAtContinuousIndex(ContinuousIndexType const & contIndex) const
{
typename TImage::IndexType index;
ConvertContinuousIndexToNearestIndex(contIndex, index);
return EvaluateAtIndex(index);
}
} // end itk namespace
#endif // NonwitnessTerminality_txx
-------------- next part --------------
A non-text attachment was scrubbed...
Name: SurfaceTerminality.h
Type: text/x-chdr
Size: 1830 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20070410/342124d1/SurfaceTerminality.h
-------------- next part --------------
#ifndef SurfaceTerminality_txx
#define SurfaceTerminality_txx
#include "SurfaceTerminality.h"
namespace itk
{
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
SurfaceTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::SurfaceTerminalityImageFunction()
{
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
SurfaceTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::Evaluate(PointType const & point) const
{
typename TImage::IndexType index;
ConvertPointToNearestIndex(point, index);
return EvaluateAtIndex(index);
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
SurfaceTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::EvaluateAtIndex(IndexType const & index) const
{
/* 9 planes, cf. Chris Pudney, p. 408, last paragraph of Section 3.1.2,
and Fig. 7:
if p, the center of the 26 neighbors,
has less than two foregroundneighbors on any of the 9 planes,
then p is the edge of a surface
*/
unsigned int const plane[9][8] = {
{1,4,7,10,16,19,22,25},
{3,4,5,12,14,21,22,23},
{9,10,11,12,14,15,16,17},
{0,4,8,9,17,18,22,26},
{2,4,6,11,15,20,22,24},
{0,1,2,12,14,24,25,26},
{6,7,8,12,14,18,19,20},
{0,3,6,10,16,20,23,26},
{2,5,8,10,16,18,21,24}
};
bool subimage[26];
TForegroundConnectivity const & fgc = TForegroundConnectivity::GetInstance();
for(int i=0; i<fgc.GetNumberOfNeighbors(); ++i)
{
itk::Offset<TForegroundConnectivity::Dimension> offset;
for(unsigned int j=0; j<TForegroundConnectivity::Dimension; ++j)
{
offset[j] = fgc.GetNeighborsPoints()[i][j];
}
subimage[i] = this->GetInputImage()->GetPixel(index+offset)!=
NumericTraits<typename TImage::PixelType>::Zero;
}
for (unsigned int pi=0; pi<9; pi++) { // loop over planes
unsigned int cnt=0;
for (unsigned int vi=0; vi<8 && cnt<2; vi++) { //control nodes of plane
unsigned int node = plane[pi][vi];
if (node>13)
node--;
if ( subimage[node] )
cnt++;
}
if (cnt < 2) // found plane with less than two foreground neighbors
return true;
}
return false;
}
template<typename TImage, typename TForegroundConnectivity,
typename TBackgroundConnectivity >
bool
SurfaceTerminalityImageFunction<TImage, TForegroundConnectivity,
TBackgroundConnectivity>
::EvaluateAtContinuousIndex(ContinuousIndexType const & contIndex) const
{
typename TImage::IndexType index;
ConvertContinuousIndexToNearestIndex(contIndex, index);
return EvaluateAtIndex(index);
}
}
#endif // SurfaceTerminality_txx
More information about the Insight-users
mailing list