[Insight-users] Edge Linking Problem
Zachary Pincus
zpincus at stanford.edu
Tue, 16 Mar 2004 09:53:02 -0800
--Apple-Mail-10-220650772
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
Hi -
I've been working on a similar problem. I've not found a perfect
solution, but here's what I've been playing with.
1) Extract Canny edges
2) Label edges -- this gets all of the connected components of a given
edge segment, but if an edge is broken into two separate segments,
there's not much to do here.
3) Find centroids of each labeled edge segment.
4) Use centroids as seed points for a level set method. Use, say,
Geodesic Active Contours, and make the feature image a function of the
Canny edge image (maybe something involving a distance map...)
With a strong curvature and edge advection term, the level set should
stick to the Canny edges, and then pretty much just bridge the gaps
otherwise.
If you wanted more precision, you could make the feature image a
function of both the original image and the Canny edges. This would
have the level set use the Canny edges when present, and the input
image when the edges are not present.
For this all to work, you'll need a ConnectedComponentImageFilter that
can label 8-connected components (as Canny edges are), and not just
face-connected (4-connected) blobs, which is what the current filter
does. I actually just finished work on such a filter. The code is
attached to this email, and as soon as I have test code worked out,
I'll submit it for inclusion in ITK.
I also found it very useful to modify the Canny filter to take as input
not simply intensity values of the edge image, but percentiles of such,
so that you don't need to re-calibrate the thresholds for each new
image. I can send that code, too, if you want.
Zach Pincus
Department of Biochemistry and Program in Biomedical Informatics
Stanford University School of Medicine
--Apple-Mail-10-220650772
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
x-mac-type=54455854;
x-unix-mode=0644;
name="itkConnectedComponentImageFilter.h"
Content-Disposition: attachment;
filename=itkConnectedComponentImageFilter.h
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkConnectedComponentImageFilter.h,v $
Language: C++
Date: $Date: 2003/09/10 14:28:45 $
Version: $Revision: 1.3 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __itkConnectedComponentImageFilter_h
#define __itkConnectedComponentImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkImage.h"
namespace itk
{
/**
* \class ConnectedComponentImageFilter
* \brief Label the objects in a binary image
*
* ConnectedComponentImageFilter labels the objects in a binary
* image. Each distinct object is assigned a unique label. The filter
* makes three passes through the image. The first pass initialized
* the output. The second pass labels each foreground pixel such that
* all the pixels associated with an object either have the same label
* or have had their labels entered into a equivalency table. The
* third pass through the image flattens the equivalency table such
* that all pixels for an object have the same label.
*
* \sa ImageToImageFilter
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT ConnectedComponentImageFilter :
public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
/**
* Standard "Self" & Superclass typedef.
*/
typedef ConnectedComponentImageFilter Self;
typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
/**
* Types from the Superclass
*/
typedef typename Superclass::InputImagePointer InputImagePointer;
/**
* Extract some information from the image types. Dimensionality
* of the two images is assumed to be the same.
*/
typedef typename TOutputImage::PixelType OutputPixelType;
typedef typename TOutputImage::InternalPixelType OutputInternalPixelType;
typedef typename TInputImage::PixelType InputPixelType;
typedef typename TInputImage::InternalPixelType InputInternalPixelType;
itkStaticConstMacro(ImageDimension, unsigned int,
TOutputImage::ImageDimension);
/**
* Image typedef support
*/
typedef TInputImage InputImageType;
typedef TOutputImage OutputImageType;
typedef typename TInputImage::IndexType IndexType;
typedef typename TInputImage::SizeType SizeType;
typedef typename TOutputImage::RegionType RegionType;
typedef std::list<IndexType> ListType;
/**
* Smart pointer typedef support
*/
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/**
* Run-time type information (and related methods)
*/
itkTypeMacro(ConnectedComponentImageFilter, ImageToImageFilter);
/**
* Method for creation through the object factory.
*/
itkNewMacro(Self);
bool GetFullyConnected() { return m_fullyConnected;}
void SetFullyConnected() { m_fullyConnected = true;}
bool GetFaceConnected() { return !m_fullyConnected;}
void SetFaceConnected() { m_fullyConnected = false;}
protected:
ConnectedComponentImageFilter() { m_fullyConnected = false; }
virtual ~ConnectedComponentImageFilter() {}
ConnectedComponentImageFilter(const Self&) {}
/**
* Standard pipeline method.
*/
void GenerateData();
/** ConnectedComponentImageFilter needs the entire input. Therefore
* it must provide an implementation GenerateInputRequestedRegion().
* \sa ProcessObject::GenerateInputRequestedRegion(). */
void GenerateInputRequestedRegion();
/** ConnectedComponentImageFilter will produce all of the output.
* Therefore it must provide an implementation of
* EnlargeOutputRequestedRegion().
* \sa ProcessObject::EnlargeOutputRequestedRegion() */
void EnlargeOutputRequestedRegion(DataObject *itkNotUsed(output));
private:
bool m_fullyConnected;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkConnectedComponentImageFilter.txx"
#endif
#endif
--Apple-Mail-10-220650772
Content-Transfer-Encoding: 7bit
Content-Type: application/text;
x-mac-type=54455854;
x-unix-mode=0644;
name="itkConnectedComponentImageFilter.txx"
Content-Disposition: attachment;
filename=itkConnectedComponentImageFilter.txx
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkConnectedComponentImageFilter.txx,v $
Language: C++
Date: $Date: 2003/09/10 14:28:45 $
Version: $Revision: 1.3 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itkConnectedComponentImageFilter_txx
#define _itkConnectedComponentImageFilter_txx
#include "itkConnectedComponentImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "itkEquivalencyTable.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkConstantBoundaryCondition.h"
namespace itk
{
template< class TInputImage, class TOutputImage >
void
ConnectedComponentImageFilter< TInputImage, TOutputImage >
::GenerateInputRequestedRegion()
{
// call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// We need all the input.
InputImagePointer input = const_cast<InputImageType *>(this->GetInput());
input->SetRequestedRegion( input->GetLargestPossibleRegion() );
}
template <class TInputImage, class TOutputImage>
void
ConnectedComponentImageFilter<TInputImage, TOutputImage>
::EnlargeOutputRequestedRegion(DataObject *)
{
this->GetOutput()
->SetRequestedRegion( this->GetOutput()->GetLargestPossibleRegion() );
}
template< class TInputImage, class TOutputImage >
void
ConnectedComponentImageFilter< TInputImage, TOutputImage >
::GenerateData()
{
// create an equivalency table
EquivalencyTable::Pointer eqTable = EquivalencyTable::New();
OutputPixelType label, originalLabel, neighborLabel;
OutputPixelType maxLabel = NumericTraits<OutputPixelType>::Zero;
const OutputPixelType maxPossibleLabel=NumericTraits<OutputPixelType>::max();
typename TOutputImage::Pointer output = this->GetOutput();
typename TInputImage::ConstPointer input = this->GetInput();
// Allocate the output and initialize to zeros
this->AllocateOutputs();
output->FillBuffer( NumericTraits<OutputPixelType>::Zero );
// Set up the boundary condition to be zero padded (used on output image)
ConstantBoundaryCondition<TOutputImage> BC;
BC.SetConstant(NumericTraits<OutputPixelType>::Zero);
// Neighborhood iterator. Let's use a shaped neighborhood so we can
// restrict the access to face connected neighbors. This iterator
// will be applied to the output image
typedef ConstShapedNeighborhoodIterator<TOutputImage> NeighborhoodIteratorType;
SizeType kernelRadius;
kernelRadius.Fill(1);
NeighborhoodIteratorType nit;
nit = NeighborhoodIteratorType(kernelRadius, output,
output->GetRequestedRegion());
nit.OverrideBoundaryCondition(&BC); // assign the boundary condition
unsigned int d;
typename NeighborhoodIteratorType::OffsetType offset;
if (m_fullyConnected)
{
// We're looking for fully-connected components:
// Activate *all* indices that are "previous" to the current
// pixel. (Previous is defined as " neighborhood index < the (0,0)
// pixel's index," which corresponds to "pixels already visited by
// the iterator.
offset.Fill(0);
int centerIndex = nit.GetNeighborhoodIndex(offset);
for (int i = 0; i < centerIndex; i++)
{
offset = nit.GetOffset(i);
nit.ActivateOffset(offset);
}
}
else
{
// OK, we're only looking for face-connected components:
// Only activate the indices that are "previous" to the current
// pixel and face connected (exclude the center pixel from the
// neighborhood)
offset.Fill(0);
for (d=0; d < InputImageType::ImageDimension; ++d)
{
offset[d] = -1;
nit.ActivateOffset(offset);
offset[d] = 0;
}
}
// along with a neighborhood iterator on the input, use a standard
// iterator on the input and output
ImageRegionConstIterator<InputImageType> it;
ImageRegionIterator<OutputImageType> oit;
it = ImageRegionConstIterator<InputImageType>(input,
output->GetRequestedRegion());
oit = ImageRegionIterator<OutputImageType>(output,
output->GetRequestedRegion());
// Setup a progress reporter. We have 3 stages to the algorithm so
// pretend we have 3 times the number of pixels
ProgressReporter progress(this, 0,
3*output->GetRequestedRegion().GetNumberOfPixels());
// Mark the output image as either background or unlabeled
it.GoToBegin();
oit.GoToBegin();
while (!it.IsAtEnd())
{
if (it.Get() != NumericTraits<InputPixelType>::Zero)
{
// mark pixel as unlabeled
oit.Set(maxPossibleLabel);
}
++it;
++oit;
progress.CompletedPixel();
}
// iterate over the image, labeling the objects and defining
// equivalence classes. Use the neighborhood iterator to access the
// "previous" neighbor pixels and an output iterator to access the
// current pixel
nit.GoToBegin();
oit.GoToBegin();
while ( !oit.IsAtEnd() )
{
// Get the current pixel label
label = oit.Get();
originalLabel = label;
// If the pixel is not background
if (label != NumericTraits<OutputPixelType>::Zero)
{
// loop over the "previous" neighbors to find labels. this loop
// may establish one or more new equivalence classes
typename NeighborhoodIteratorType::ConstIterator sIt;
for (sIt = nit.Begin(); !sIt.IsAtEnd(); ++sIt)
{
// get the label of the pixel previous to this one along a
// particular dimension (neighbors activated in neighborhood iterator)
neighborLabel = sIt.Get();
// if the previous pixel has a label, verify equivalence or
// establish a new equivalence
if (neighborLabel != NumericTraits<OutputPixelType>::Zero)
{
// if current pixel is unlabeled, copy the label from neighbor
if (label == maxPossibleLabel)
{
// copy the label from the previous pixel
label = neighborLabel;
}
// else if current pixel has a label that is not already
// equivalent to the label of the previous pixel, then setup
// a new equivalence. note the use of Lookup() and not
// RecursiveLookup(). this is possible since we keep the
// equivalence table flat.
else if ((label != neighborLabel)
&& (eqTable->Lookup(label) != eqTable->Lookup(neighborLabel)))
{
eqTable->Add(label, neighborLabel);
// if we keep the equivalency table up to date, then we
// can use straight calls to Lookup() instead of
// RecursiveLookUp(). This works out to be 3X faster.
eqTable->Flatten();
}
}
}
// if none of the "previous" neighbors were set, then make a new label
if (originalLabel == label)
{
// create a new entry label
if (maxLabel == maxPossibleLabel)
{
itkWarningMacro(<< "ConnectedComponentImageFilter::GenerateData: Number of labels exceeds number of available labels for the output type." );
}
else
{
++maxLabel;
}
// assign the new label
label = maxLabel;
}
// Finally, set the output pixel to whatever label we have
if (label != originalLabel)
{
oit.Set( label );
}
}
// move the iterators
++nit;
++oit;
progress.CompletedPixel();
}
// Flatten the equavalency table
eqTable->Flatten();
// remap the labels
oit.GoToBegin();
while ( !oit.IsAtEnd() )
{
label = oit.Get();
// if pixel has a label, write out the final equivalence
if (label != NumericTraits<OutputPixelType>::Zero)
{
oit.Set( eqTable->Lookup( label ) );
}
++oit;
progress.CompletedPixel();
}
}
} // end namespace itk
#endif
--Apple-Mail-10-220650772
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain;
charset=ISO-8859-1;
format=flowed
On Mar 12, 2004, at 10:19 PM, CSPL wrote:
> Hello Group,
> =A0
> =A0I have a requirement of segmenting(extracting) an object of =
interest=20
> from an image and the method I am supposed to use is the following:
> =A0
> =A01. Find the edges of the required object.
> =A02. Given the edges of the object, fill or some how form the =
complete=20
> object.
> =A0=A0i.e, define all pixels of the object whose edges we have found =
in=20
> step 1.
> =A0
> =A0To achieve this, I used Canny Edge Detection filter for finding the=20=
> edges.
> =A0Here is the actual problem:
> =A0The edges formed by the filter are not connected. The object is=20
> surrounded by broken=A0 edge lines.
> =A0
> =A0How do I get the completely connected edge of the object? If there =
is=20
> a way to link these broken edges, on what factor will the accuracy of=20=
> this edge linking depend?
> =A0
> =A0The documentation for Canny Edge Detection Filter says the =
following:
> =A0"TODO : Edge-linking will be added when an itk connected component=20=
> labeling algorithm is available."
> =A0
> =A0"ConnectedComponentImageFilter" is already available in ITK. Is =
this=20
> the class being referred to in the above line?=A0If it is the same, =
then=20
> how can I use it in this situation?
> =A0
> =A0Is it possible to move to the second step without obtaining the=20
> complete edge of the object?
> =A0
> Regards,
> Rama Krishna.
--Apple-Mail-10-220650772--