Summary:
The authors discuss how to extract regional minima/maxima regions from binary/greyscale images. Two approaches are presented and compared: 1. the current ITK approach using the HConvexImageFilter and HConcaveImageFilter, and 2. a new approach using a simple flooding-based algorithm (referred to by the authors as ValuedRegionalExtremaImageFilter).
Hypothesis:
The hypothesis is that the ValuedRegionalExtremaImageFilter approach is computational faster than the HConvexImageFilter/HConcaveImageFilter approach. Test results indicate that the proposed approach (ValuedRegionalExtremaImageFilter) is faster.
Evidence:
Each approach has been implemented in ITK and execution times reported in the paper.
Open Science:
Full open source for the filters, examples and tests, as well as input images, are provided.
Reproducibility:
I was able to download, compile, and run the work. However, I had some minor problems compiling the code because I was using an older version of ITK (InsightToolkit-2.2.0). Compiling with InsightToolkit-2.4.1 was successful. I was also unable to run the tests because they rely on the DART Image Compare program. (Does anyone know how I can obtain the Image Compare program for testing Insight Journal articles? The IJ WIKI mentions that DART tests use Image Compare, but it does not explain how to obtain it for desktop testing/reproducibility.)
Use of Open Source Software:
Uses ITK.
Open Source Contributions:
The source code was relatively easy to get going (besides the minor unrelated issues listed above).
Code Quality:
The code quality is fairly easy to read. However, I would like to see more comments within ValuedRegionalExtremaImageFilter::GenerateData() more fully explaining the algorithm - I find it hard to read what is being done and why.
Final Thoughts:
Thanks for this addition to ITK. The speed increases for finding regional maxima/minima will help enhance various segmentation tasks. A comprehensive, well presented paper.
Comment by Dan Mueller: Additional comments, and flood filling
I like the addition comments you have added so far - I found it much easier to work out how the algorithm works.
Now that I better understand what is happening, have you considered using the itk::FloodFilledFunctionConditionalIterator? This iterator is templated over a comparison function and allows you to walk over a neighbourhood image region which conforms to the comparison function. The only drawback is this iterator is currently *not* implemented with a ShapedNeighborhoodIterator and therefore you can not set the FullyConnected property... From what I understand however, it does what you are doing in the last phase of your algorithm... Just a thought.
Comment by Richard Beare: Draft of additional comments
Hi Dan,
I've started extending the comments in ValuedRegionalExtremaImageFilter::GenerateData() - see below. Is there any particular part where you think more detail is helpful?
ValuedRegionalExtremaImageFilter
::GenerateData()
{
// Allocate the output
this->AllocateOutputs();
// 2 phases
ProgressReporter progress(this, 0, this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()*2);
// copy input to output - isn't there a better way?
typedef ImageRegionConstIterator InputIterator;
typedef ImageRegionIterator OutputIterator;
InputIterator inIt( this->GetInput(),
this->GetOutput()->GetRequestedRegion() );
OutputIterator outIt( this->GetOutput(),
this->GetOutput()->GetRequestedRegion() );
inIt = inIt.Begin();
outIt = outIt.Begin();
InputImagePixelType firstValue = inIt.Get();
this->m_Flat = true;
while ( !outIt.IsAtEnd() )
{
InputImagePixelType currentValue = inIt.Get();
outIt.Set( static_cast( currentValue ) );
if( currentValue != firstValue )
{ this->m_Flat = false; }
++inIt;
++outIt;
progress.CompletedPixel();
}
// if the image is flat, there is no need to do the work:
// the image will be unchanged
if( !this->m_Flat )
{
// Now for the real work!
// More iterators - use shaped ones so that we can set connectivity
// Note : all comments refer to finding regional minima, because
// it is briefer and clearer than trying to describe both regional
// maxima and minima processes at the same time
ISizeType kernelRadius;
kernelRadius.Fill(1);
NOutputIterator outNIt(kernelRadius,
this->GetOutput(),
this->GetOutput()->GetRequestedRegion() );
setConnectivity( &outNIt, m_FullyConnected );
CNInputIterator inNIt(kernelRadius,
this->GetInput(),
this->GetOutput()->GetRequestedRegion() );
setConnectivity( &inNIt, m_FullyConnected );
ConstantBoundaryCondition iBC;
iBC.SetConstant(m_MarkerValue);
inNIt.OverrideBoundaryCondition(&iBC);
ConstantBoundaryCondition oBC;
oBC.SetConstant(m_MarkerValue);
outNIt.OverrideBoundaryCondition(&oBC);
TFunction1 compareIn;
TFunction2 compareOut;
outIt = outIt.Begin();
// set up the stack and neighbor list
IndexStack IS;
typename NOutputIterator::IndexListType IndexList;
IndexList = outNIt.GetActiveIndexList();
while ( !outIt.IsAtEnd() )
{
OutputImagePixelType V = outIt.Get();
// if the output pixel value = the marker value then we have
// already visited this pixel and don't need to do so again
if (compareOut(V, m_MarkerValue))
{
// reposition the input iterator
inNIt += outIt.GetIndex() - inNIt.GetIndex();
InputImagePixelType Cent = static_cast(V);
// check each neighbor of the input pixel
typename CNInputIterator::ConstIterator sIt;
for (sIt = inNIt.Begin(); !sIt.IsAtEnd(); ++sIt)
{
InputImagePixelType Adjacent = sIt.Get();
if (compareIn(Adjacent, Cent))
{
// The centre pixel cannot be part of a regional minima
// because one of its neighbors is smaller.
// Set all pixels in the output image that are connected to
// the centre pixel and have the same value to
// m_MarkerValue
// This is the flood filling step. It is a simple, stack
// based, procedure. The original value (V) of the pixel is
// recorded and the pixel index in the output image
// is set to the marker value. The stack is initialized
// with the pixel index. The flooding procedure pops the
// stack, sets that index to the marker value and places the
// indexes of all neighbors with value V on the stack. The
// process terminates when the stack is empty.
outNIt += outIt.GetIndex() - outNIt.GetIndex();
//setConnectedPixels(outNIt, V, IS, IndexList);
OutputImagePixelType NVal;
//IndexStack IS;
OutIndexType idx;
// Initialize the stack
IS.push(outNIt.GetIndex());
outNIt.SetCenterPixel(m_MarkerValue);
typename NOutputIterator::IndexListType::const_iterator LIt;
while (!IS.empty())
{
// Pop the stack
idx = IS.top();
IS.pop();
// position the iterator
outNIt += idx - outNIt.GetIndex();
// check neighbors
for (LIt = IndexList.begin(); LIt != IndexList.end(); ++LIt)
{
NVal = outNIt.GetPixel(*LIt);
if (NVal == V)
{
// still in a flat zone
IS.push(outNIt.GetIndex(*LIt));
// set the output as the marker value
outNIt.SetPixel(*LIt, m_MarkerValue);
}
}
}
// end flooding
}
}
}
++outIt;
progress.CompletedPixel();
}
}
}
Comment by Gaetan Lehmann: ImageCompare
Hi Dan,
You can get ImageCompare at http://insight-journal.org/documentation/CMake/ImageCompare.tgz
Thanks for your review.
Gaetan
2 Files (2Mb)
Automatic Testing Results
by Insight-Journal Dashboard
on Thu Jan 26 16:50:51 2006 for revision #3 
GeoInterp: Contour Interpolation with Geodesic Snakes
Geodesic Snakes
N4ITK: Nick's N3 ITK Implementation For MRI Bias Field Correction
