itkWatershedSegmenter.txx

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkWatershedSegmenter.txx,v $
  Language:  C++
  Date:      $Date: 2009-01-27 19:30:19 $
  Version:   $Revision: 1.37 $

  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 __itkWatershedSegmenter_txx
#define __itkWatershedSegmenter_txx

#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodIterator.h"
#include <stack>
#include <list>

namespace itk
{
namespace watershed
{

/*
  ----------------------------------------------------------------------------
  Algorithm methods
  ----------------------------------------------------------------------------
*/


template <class TInputImage>
Segmenter<TInputImage>::~Segmenter()
{
  if (m_Connectivity.index != 0) 
    {
    delete[] m_Connectivity.index;
    }
  if (m_Connectivity.direction !=0 ) 
    {
    delete[] m_Connectivity.direction;
    }
}


template <class TInputImage>
void Segmenter<TInputImage>::GenerateData()
{
  //
  // Allocate all the necessary temporary data structures and variables that
  // will be used in this algorithm.  Also re-initialize some temporary data
  // structures that may have been used in previous updates of this filter.
  //
  unsigned int i;

  this->UpdateProgress(0.0);
  if (m_DoBoundaryAnalysis == false)
    {
    this->GetSegmentTable()->Clear();
    this->SetCurrentLabel(1);
    }
  
  flat_region_table_t flatRegions;

  typename InputImageType::Pointer input   = this->GetInputImage();
  typename OutputImageType::Pointer output = this->GetOutputImage();
  typename BoundaryType::Pointer boundary  = this->GetBoundary();

  // ------------------------------------------------------------------------
  //
  // HERE ARE THE ASSUMPTIONS ABOUT REGION SIZES FOR NOW.  WHEN THE PIPELINE
  // FULLY SUPPORTS STREAMING, THESE WILL NEED TO BE CHANGED ACCORDINGLY.
  //
  // 1) All region sizes are equivalent.  There is no distinction among
  // regions.  The region size is assumed to be padded one pixel out along each 
  // chunk face unless that face touches an actual data set boundary.
  //
  // 2) The ivar m_LargestPossibleRegion represents the actual size of the data
  // set.  This has to be set by the user since the pipeline sometimes clobbers
  // the actual LargestPossibleRegion (?).
  //
  // -------------------------------------------------------------------------
  
  //
  // Generate the "face" regions A that constitute our shared boundary with
  // another chunk.  Also determine which face regions B lie on a the true
  // dataset boundary.  The faces corresponding to B will need to be padded
  // out a pixel when we threshold so that we can construct the retaining wall
  // along those faces.
  //
  ImageRegionType regionToProcess       = output->GetRequestedRegion();
  ImageRegionType largestPossibleRegion = this->GetLargestPossibleRegion();
  ImageRegionType thresholdImageRegion  = regionToProcess;
  ImageRegionType thresholdLargestPossibleRegion
    = this->GetLargestPossibleRegion();
  
  // First we have to find the boundaries and adjust the threshold image size
  typename ImageRegionType::IndexType  tidx= thresholdImageRegion.GetIndex();
  typename ImageRegionType::SizeType   tsz = thresholdImageRegion.GetSize();
  typename ImageRegionType::IndexType  tlidx= thresholdLargestPossibleRegion.GetIndex();
  typename ImageRegionType::SizeType   tlsz = thresholdLargestPossibleRegion.GetSize();
  for (i = 0; i < ImageDimension; ++i)
    {
    ImageRegionType reg;
    typename ImageRegionType::IndexType idx = regionToProcess.GetIndex();
    typename ImageRegionType::SizeType  sz  = regionToProcess.GetSize();

    // Set LOW face
    idx[i] = regionToProcess.GetIndex()[i];
    sz[i]  = 1;
    reg.SetSize(sz);
    reg.SetIndex(idx);

    if (reg.GetIndex()[i] == largestPossibleRegion.GetIndex()[i])
      {
      // This is facing a true data set boundary
      tsz[i] += 1; // we need to pad our threshold image on this face
      tidx[i] -= 1;
      tlsz[i] += 1; // we need to pad our threshold image on this face
      tlidx[i] -= 1;
      
      boundary->SetValid(false, i, 0);
      }
    else
      {
      // This is an overlap with another data chunk in the data set
      // Mark this boundary face as valid.
      boundary->SetValid(true, i, 0);
      }

    // Set HIGH face
    idx[i] = (regionToProcess.GetIndex()[i]+regionToProcess.GetSize()[i]) - 1;
    reg.SetSize(sz);
    reg.SetIndex(idx);
    if ( (reg.GetIndex()[i] + reg.GetSize()[i])
         == (largestPossibleRegion.GetIndex()[i]
             + largestPossibleRegion.GetSize()[i]) )
      {
      // This is facing a true data set boundary
      tsz[i] += 1; // we need to pad our threshold image on this face
      tlsz[i] += 1; // we need to pad our threshold image on this face
      boundary->SetValid(false, i, 1);
      }
    else
      {
      // This is an overlap with another data chunk in the data set
      // Mark this face as valid in the boundary.
      boundary->SetValid(true, i, 1);
      }
    }
  thresholdImageRegion.SetSize(tsz);
  thresholdImageRegion.SetIndex(tidx);
  thresholdLargestPossibleRegion.SetSize(tlsz);
  thresholdLargestPossibleRegion.SetIndex(tlidx);

  // Now create and allocate the threshold image.  We need a single pixel
  // border around the NxM region we are segmenting.  This means that for faces 
  // that have no overlap into another chunk, we have to pad the image.
  typename InputImageType::Pointer thresholdImage = InputImageType::New();

  thresholdImage->SetLargestPossibleRegion(thresholdLargestPossibleRegion);
  thresholdImage->SetBufferedRegion(thresholdImageRegion);
  thresholdImage->SetRequestedRegion(thresholdImageRegion);
  thresholdImage->Allocate();

  // Now threshold the image. First we calculate the dynamic range of
  // the input.  Then, the threshold operation clamps the lower
  // intensity values at the prescribed threshold.  If the data is
  // integral, then any intensity at NumericTraits<>::max() is reduced
  // by one intensity value.  This allows the watershed algorithm to
  // build a barrier around the image with values above the maximum
  // intensity value which trivially stop the steepest descent search
  // for local minima without requiring expensive boundary conditions.
  //
  //
  InputPixelType minimum, maximum;
  Self::MinMax(input, regionToProcess, minimum, maximum);
  // cap the maximum in the image so that we can always define a pixel
  // value that is one greater than the maximum value in the image.
  if (NumericTraits<InputPixelType>::is_integer
      && maximum == NumericTraits<InputPixelType>::max())
    {
    maximum -= NumericTraits<InputPixelType>::One;
    }
  // threshold the image. 
  Self::Threshold(thresholdImage, input, regionToProcess, regionToProcess,
                  static_cast<InputPixelType>((m_Threshold * (maximum - minimum)) + minimum));

  
  //
  // Redefine the regionToProcess in terms of the threshold image.  The region
  // to  process represents all the pixels contained within the 1 pixel padded
  // boundary of the threshold image.
  //
  typename ImageRegionType::SizeType  irsz;
  typename ImageRegionType::IndexType iridx;
  for (i = 0; i < ImageDimension; ++i)
    {
    irsz[i]  = thresholdImageRegion.GetSize()[i] - 2;
    iridx[i] = thresholdImageRegion.GetIndex()[i] + 1;
    }
  regionToProcess.SetIndex(iridx);
  regionToProcess.SetSize(irsz);

  //
  // Initialize the connectivity information that will be used by the 
  // segmentation algorithm.
  //
  this->GenerateConnectivity();

  //
  // Store the regionToProcess in the RequestedRegion of the threshold image.
  // We are now completely done with the input image.  The input image memory
  // can be released at this point if need be.
  //
  thresholdImage->SetRequestedRegion(regionToProcess);
  this->ReleaseInputs();
  
  //
  // At this point we are ready to define the output
  // buffer and allocate memory for the output image.
  //
  output->SetBufferedRegion(thresholdImage->GetBufferedRegion());
  output->Allocate();
  Self::SetOutputImageValues(output, output->GetBufferedRegion(), Self::NULL_LABEL);

  //
  // Now we can create appropriate boundary regions for analyzing the
  // flow at the boundaries from the requested region of the threshold
  // image.
  //
  typename BoundaryType::IndexType b_idx;
  ImageRegionType reg_b;
  typename ImageRegionType::IndexType idx_b;
  typename ImageRegionType::SizeType  sz_b;

  for (b_idx.first = 0; b_idx.first < ImageDimension; ++b_idx.first)
    {
    for (b_idx.second = 0; b_idx.second < 2; ++b_idx.second)
      {
      if (boundary->GetValid(b_idx) == false) continue;
      idx_b = thresholdImage->GetRequestedRegion().GetIndex();
      sz_b = thresholdImage->GetRequestedRegion().GetSize();

      if (b_idx.second == 1) // HIGH face must adjust start index
        idx_b[b_idx.first] += sz_b[b_idx.first] - 1;

      sz_b[b_idx.first] = 1;

      reg_b.SetIndex(idx_b);
      reg_b.SetSize(sz_b);
        
      boundary->GetFace(b_idx)->SetLargestPossibleRegion(reg_b);
      boundary->GetFace(b_idx)->SetRequestedRegion(reg_b);
      boundary->GetFace(b_idx)->SetBufferedRegion(reg_b);
      boundary->GetFace(b_idx)->Allocate();
      }
    }
  this->UpdateProgress(0.1);
  
  
  //
  // Analyze the flow at the boundaries.  This method labels all the boundary
  // pixels that flow out of this chunk (either through gradient descent or
  // flat-region connectivity) and constructs the appropriate Boundary 
  // data structures.
  //
  if (m_DoBoundaryAnalysis == true)
    {
    this->InitializeBoundary();
    this->AnalyzeBoundaryFlow(thresholdImage, flatRegions, maximum +
                              NumericTraits<InputPixelType>::One);
    }

  this->UpdateProgress(0.2);

  //
  // Build a ``retaining wall'' around the image so that gradient descent
  // analysis can be done without worrying about boundaries.
  // 
  // All overlap boundary information will be overwritten, but is no longer
  // needed now.
  // 
  this->BuildRetainingWall( thresholdImage,
                            thresholdImage->GetBufferedRegion(),
                            maximum + NumericTraits<InputPixelType>::One );

  //
  // Label all the local minima pixels in the image.  This function also
  // labels flat regions, defined as regions where connected pixels all have
  // the same value.
  //
  this->LabelMinima(thresholdImage, thresholdImage->GetRequestedRegion(),
                    flatRegions, maximum + NumericTraits<InputPixelType>::One);
  this->UpdateProgress(0.3);

  this->GradientDescent(thresholdImage, thresholdImage->GetRequestedRegion());
  this->UpdateProgress(0.4);

  this->DescendFlatRegions(flatRegions, thresholdImage->GetRequestedRegion());
  this->UpdateProgress(0.5);

  this->UpdateSegmentTable(thresholdImage, thresholdImage->GetRequestedRegion());
  this->UpdateProgress(0.6);

  if (m_DoBoundaryAnalysis == true)
    {  this->CollectBoundaryInformation(flatRegions); }
  this->UpdateProgress(0.7);

  if (m_SortEdgeLists == true)
    {  this->GetSegmentTable()->SortEdgeLists(); }
  this->UpdateProgress(0.8);

  this->GetSegmentTable()->SetMaximumDepth(maximum - minimum);
  this->UpdateProgress(1.0);

}

template <class TInputImage>
void Segmenter<TInputImage>  
::CollectBoundaryInformation(flat_region_table_t &flatRegions)
{
  typename OutputImageType::Pointer output = this->GetOutputImage();
  typename BoundaryType::Pointer boundary  = this->GetBoundary();

  ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t> faceIt;
  ImageRegionIterator<OutputImageType>              labelIt;

  typename BoundaryType::face_t::Pointer        face;
  typedef typename BoundaryType::flat_hash_t    flats_t;
  typename BoundaryType::flat_hash_t           *flats;
  typename BoundaryType::flat_hash_t::iterator  flats_it;
  typename BoundaryType::flat_region_t          flr;
  typename flat_region_table_t::iterator        flrt_it;
  
  typename BoundaryType::IndexType       idx;
  ImageRegionType                        region;

  for (idx.first = 0; idx.first < ImageDimension; (idx.first)++)
    {
    for (idx.second = 0; idx.second < 2; (idx.second)++)
      {
      if (boundary->GetValid(idx) == false) continue;
        
      face  = boundary->GetFace(idx);
      flats = boundary->GetFlatHash(idx);
      region = face->GetRequestedRegion();
        
      // Grab all the labels of the boundary pixels.
      faceIt = ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t> (face,
                                                                       region);
      labelIt = ImageRegionIterator<OutputImageType> (output, region);
      faceIt = faceIt.Begin();
      labelIt = labelIt.Begin();
      while (! faceIt.IsAtEnd() )
        {
        faceIt.Value().label = labelIt.Get();
            
        // Is this a flat region that flows out?
        flrt_it = flatRegions.find(labelIt.Get());
        if ( faceIt.Get().flow != NULL_FLOW
             && flrt_it != flatRegions.end() ) 
          {
          // Have we already entered this
          // flat region into the boundary?
          flats_it = flats->find(labelIt.Get());
          if (flats_it == flats->end())  // NO
            {
            flr.bounds_min = (*flrt_it).second.bounds_min;
            flr.min_label  = *((*flrt_it).second.min_label_ptr);
            flr.value      = (*flrt_it).second.value;
            flr.offset_list.push_back(
              face->ComputeOffset(faceIt.GetIndex()));
            flats->insert(
              BoundaryFlatHashValueType(labelIt.Get(), flr));
            flr.offset_list.clear();
            }
          else // YES
            {
            (*flats_it).second.offset_list.push_back(face->ComputeOffset(faceIt.GetIndex()));
            }
          }
            
        ++faceIt;
        ++labelIt;
        }
      }
    }
}
  
template <class TInputImage>
void Segmenter<TInputImage>
::InitializeBoundary()
{
  ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t> faceIt;
  typename BoundaryType::face_t::Pointer face;
  typename BoundaryType::face_pixel_t fps;
  BoundaryIndexType idx;
  
  fps.flow  = NULL_FLOW;
  fps.label = NULL_LABEL;
  
  for (idx.first = 0; idx.first < ImageDimension; ++(idx.first))
    {
    for (idx.second = 0; idx.second < 2; ++(idx.second))
      {
      if (this->GetBoundary()->GetValid(idx) == false) continue;
      this->GetBoundary()->GetFlatHash(idx)->clear();
      face = this->GetBoundary()->GetFace(idx);
      faceIt = ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t>
        (face, face->GetBufferedRegion());
      for (faceIt = faceIt.Begin(); ! faceIt.IsAtEnd(); ++faceIt)
        faceIt.Set(fps);
      }
    }
}
  
template <class TInputImage>
void Segmenter<TInputImage>
::AnalyzeBoundaryFlow(InputImageTypePointer thresholdImage,
                      flat_region_table_t &flatRegions,
                      InputPixelType max)
{
  //
  // NOTE: For ease of initial implementation, this method does
  // not support arbitrary connectivity across boundaries (yet). 10-8-01 jc
  //
  unsigned int nCenter, i, nPos, cPos;
  bool isSteepest;
  ConstNeighborhoodIterator<InputImageType>          searchIt;
  NeighborhoodIterator<OutputImageType>               labelIt;
  ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t> faceIt;

  BoundaryIndexType idx;
  ImageRegionType region;
  typename ConstNeighborhoodIterator<InputImageType>::RadiusType rad;

  typename BoundaryType::face_pixel_t   fps;
  flat_region_t tempFlatRegion;

  typename OutputImageType::Pointer output   = this->GetOutputImage();
  typename BoundaryType::Pointer    boundary = this->GetBoundary();
  
  for (i = 0; i < ImageDimension; ++i)
    {
    rad[i] = 1;
    }
  fps.label = NULL_LABEL;

  EquivalencyTable::Pointer eqTable = EquivalencyTable::New();

  
  // Process each boundary region.
  for (idx.first = 0; idx.first < ImageDimension; ++(idx.first))
    {
    for (idx.second = 0; idx.second < 2; ++(idx.second))
      {
      // Skip irrelevant boundaries
      if (boundary->GetValid(idx) == false) continue;
          
      typename BoundaryType::face_t::Pointer face = boundary->GetFace(idx);
      region = face->GetRequestedRegion();

      searchIt
        = ConstNeighborhoodIterator<InputImageType> (rad, thresholdImage, region);
      labelIt = NeighborhoodIterator<OutputImageType> (rad, output, region);
      faceIt  = ImageRegionIterator<ITK_TYPENAME BoundaryType::face_t> (face, region);
          
      nCenter = searchIt.Size() / 2;
      searchIt.GoToBegin();
      labelIt.GoToBegin();
          
      if ((idx).second == 0)
        {
        // Low face
        cPos = m_Connectivity.index[(idx).first];
        }
      else
        {
        // High face
        cPos = m_Connectivity.index[(ImageDimension - 1)
                                    + (ImageDimension - (idx).first)];
        }

      while ( ! searchIt.IsAtEnd() )
        {
        // Is this a flat connection?
        if ( searchIt.GetPixel(nCenter) == searchIt.GetPixel(cPos) )
          {
          // Fill in the boundary flow information.
          // Labels will be collected later.
          fps.flow   = static_cast<short>(cPos);
          faceIt.Set(fps);

          // Are we touching flat regions
          // that have already been labeled?
          bool _labeled    = false;
          bool _connected  = false;
          for ( i=0; i <m_Connectivity.size; i++)
            {
            nPos = m_Connectivity.index[i];
            if (   searchIt.GetPixel(nCenter) == searchIt.GetPixel(nPos)
                   && labelIt.GetPixel(nPos) != Self::NULL_LABEL
                   && labelIt.GetPixel(nPos) != labelIt.GetPixel(nCenter)
              )
              {
              _connected = true;
              if (_labeled == false)
                {
                labelIt.SetPixel(nCenter,
                                 labelIt.GetPixel(nPos));
                _labeled = true;
                }
              else
                {
                eqTable->Add(labelIt.GetPixel(nCenter), labelIt.GetPixel(nPos));
                }
              }
            }
          if (_connected == false ) // Add a new flat region.
            {
            labelIt.SetPixel(nCenter, m_CurrentLabel);
                      
            // Add a flat region to the (global) flat region table
            tempFlatRegion.bounds_min    = max;
            tempFlatRegion.min_label_ptr = output->GetBufferPointer() +
              output->ComputeOffset(labelIt.GetIndex());
            tempFlatRegion.value         = searchIt.GetPixel(nCenter);
            tempFlatRegion.is_on_boundary= true;
            flatRegions[m_CurrentLabel]  = tempFlatRegion;
                      
            m_CurrentLabel++;
            }
          }
        else  // Is cPos the path of steepest descent?
          {
          if ( searchIt.GetPixel(cPos) < searchIt.GetPixel(nCenter) )
            {
            isSteepest = true;
            for (i = 0; i < m_Connectivity.size; i++)
              {
              nPos = m_Connectivity.index[i];
              if (searchIt.GetPixel(nPos) < searchIt.GetPixel(cPos))
                {
                isSteepest = false;
                break;
                }
              }
            }
          else isSteepest = false;
                  
          if (isSteepest == true)
            {
            // Label this pixel. It will be safely treated as a local
            // minimum by the rest of the segmentation algorithm.
            labelIt.SetPixel(nCenter, m_CurrentLabel);

            // Add the connectivity information
            // to the boundary data structure.
            fps.flow  = static_cast<short>(cPos);
            faceIt.Set(fps);

            // Since we've labeled this pixel, we need to check to
            // make sure this is not also a flat region.  If it is,
            // then it must be entered into the flat region table
            // or we could have problems later on.
            for (i = 0; i < m_Connectivity.size; i++)
              {
              nPos = m_Connectivity.index[i];
              if ( searchIt.GetPixel(nPos) ==
                   searchIt.GetPixel(nCenter) )
                {
                tempFlatRegion.bounds_min = max;
                tempFlatRegion.min_label_ptr =
                  output->GetBufferPointer() +
                  output->ComputeOffset(labelIt.GetIndex());
                tempFlatRegion.value =
                  searchIt.GetPixel(nCenter);
                tempFlatRegion.is_on_boundary = false;
                flatRegions[m_CurrentLabel] = tempFlatRegion;
                break;
                }
              }
            m_CurrentLabel++;
            }
          }
              
        ++searchIt;
        ++labelIt;
        ++faceIt;
        }
      }
    }

  eqTable->Flatten();

  // Now relabel any equivalent regions in the boundaries.
  for (idx.first = 0; idx.first < ImageDimension; ++(idx.first))
    {
    for (idx.second = 0; idx.second < 2; ++(idx.second))
      {
      // Skip irrelevant boundaries
      if (boundary->GetValid(idx) == false) continue;
          
      typename BoundaryType::face_t::Pointer face = boundary->GetFace(idx);
      region = face->GetRequestedRegion();

      Self::RelabelImage(output, region, eqTable);
      }
    }

  // Merge the flat regions in the table
  Self::MergeFlatRegions(flatRegions, eqTable);
}
  
template <class TInputImage>
void Segmenter<TInputImage>
::GenerateConnectivity()
{
  unsigned int i, j, nSize, nCenter, stride;
  int d;

  //
  // Creates city-block style connectivity.  4-Neighbors in 2D.  6-Neighbors in 
  // 3D, etc...  Order of creation MUST be lowest index to highest index in the 
  // neighborhood.  I.e. for 4 connectivity,
  //
  //     * 1 *
  //     2 * 3
  //     * 4 * 
  //
  // Algorithms assume this order to the connectivity.
  //
  typename ConstNeighborhoodIterator<InputImageType>::RadiusType rad;
  for (i = 0; i < ImageDimension; ++i)
    {
    rad[i] = 1;
    }
  ConstNeighborhoodIterator<InputImageType> it(rad, this->GetInputImage(),
                                               this->GetInputImage()->GetRequestedRegion());
  nSize   = it.Size();
  nCenter = nSize >> 1;
  
 
  for (i =0; i < m_Connectivity.size; i++) // initialize move list
    {
    for (j =0; j < ImageDimension; j++)
      {
      m_Connectivity.direction[i][j] = 0;
      }
    }
  i = 0;
  for (d = ImageDimension-1; d >=0; d--)
    {
    stride = it.GetStride(d);
    m_Connectivity.index[i] = nCenter - stride;
    m_Connectivity.direction[i][d] = -1;
    i++;
    }
  for (d = 0; d < static_cast<int>(ImageDimension); d++)
    {
    stride = it.GetStride(d);
    m_Connectivity.index[i] = nCenter + stride;
    m_Connectivity.direction[i][d] = 1;
    i++;
    }
}
  
template <class TInputImage>
void Segmenter<TInputImage>
::LabelMinima(InputImageTypePointer img, ImageRegionType region,
              typename Self::flat_region_table_t &flatRegions, InputPixelType Max)
{
  unsigned int i, nSize, nCenter, nPos = 0;
  bool foundSinglePixelMinimum, foundFlatRegion;
  InputPixelType maxValue = Max;
  
  flat_region_t tempFlatRegion;
  typename flat_region_table_t::iterator flatPtr;
  InputPixelType currentValue;
  EquivalencyTable::Pointer equivalentLabels = EquivalencyTable::New();

  typename OutputImageType::Pointer output = this->GetOutputImage();


  // Set up the iterators.
  typename ConstNeighborhoodIterator<InputImageType>::RadiusType rad;
  for (i = 0; i < ImageDimension; ++i)
    {
    rad[i] = 1;
    }
  ConstNeighborhoodIterator<InputImageType>   searchIt(rad, img, region);
  NeighborhoodIterator<OutputImageType>       labelIt(rad, output, region);
  nSize   = searchIt.Size();
  nCenter = nSize >> 1;

  // Sweep through the images.  Label all local minima
  // and record information for all the flat regions.
  for (searchIt.GoToBegin(), labelIt.GoToBegin();
       ! searchIt.IsAtEnd(); ++searchIt, ++labelIt)
    {
    foundSinglePixelMinimum = true;
    foundFlatRegion = false;

    // If this pixel has been labeled already,
    // skip directly to the next iteration.
    if ( labelIt.GetPixel(nCenter) != Self::NULL_LABEL )  continue;
      
    // Compare current pixel value with its neighbors.
    currentValue = searchIt.GetPixel(nCenter);

    for (i = 0; i < m_Connectivity.size; ++i)
      {
      nPos = m_Connectivity.index[i];
      if ( currentValue == searchIt.GetPixel(nPos) )
        {
        foundFlatRegion  = true;
        break;
        }
      else if ( currentValue > searchIt.GetPixel(nPos) )
        {
        foundSinglePixelMinimum = false;
        }
      }

    if (foundFlatRegion)
      {
      if ( labelIt.GetPixel(nPos) != Self::NULL_LABEL )// If the flat region is already
        {                                              // labeled, label this to match.
        labelIt.SetPixel(nCenter, labelIt.GetPixel(nPos));
        }
      else // Add a new flat region to the table.
        {  // Initialize its contents.
        labelIt.SetPixel(nCenter,  m_CurrentLabel);
        nPos = m_Connectivity.index[0];

        tempFlatRegion.bounds_min        = maxValue;
        tempFlatRegion.min_label_ptr     = labelIt[nPos];
        tempFlatRegion.value             = currentValue;
        flatRegions[m_CurrentLabel]      = tempFlatRegion;
        m_CurrentLabel = m_CurrentLabel + 1;
        }
          
      // While we're at it, check to see if we have just linked two flat
      // regions with the same height value.  Save that info for later.
      for ( i++; i <m_Connectivity.size; ++i)
        {
        nPos = m_Connectivity.index[i];
        if (   searchIt.GetPixel(nCenter) == searchIt.GetPixel(nPos)
               && labelIt.GetPixel(nPos) != Self::NULL_LABEL
               && labelIt.GetPixel(nPos) != labelIt.GetPixel(nCenter)
          )
          {
          equivalentLabels->Add(labelIt.GetPixel(nCenter),
                                labelIt.GetPixel(nPos));
          }
        }
      }
    else if (foundSinglePixelMinimum)
      {
      labelIt.SetPixel(nCenter,  m_CurrentLabel);
      m_CurrentLabel = m_CurrentLabel + 1;
      } 
    }

  // Merge the flat regions that we identified as connected components.
  Self::MergeFlatRegions(flatRegions, equivalentLabels);

  // Relabel the image with the merged regions.
  Self::RelabelImage(output, region, equivalentLabels);

  equivalentLabels->Clear();
  
  // Now make another pass to establish the
  // boundary values for the flat regions.
  for (searchIt.GoToBegin(), labelIt.GoToBegin();
       ! searchIt.IsAtEnd(); ++searchIt, ++labelIt)
    {
    flatPtr = flatRegions.find( labelIt.GetPixel(nCenter) ); 
    if (  flatPtr != flatRegions.end() ) // If we are in a flat region
      { // Search the connectivity neighborhood
      // for lesser boundary pixels.
      for (i = 0; i < m_Connectivity.size; ++i)
        {
        nPos = m_Connectivity.index[i];

        if  ( labelIt.GetPixel(nPos) != labelIt.GetPixel(nCenter) &&
              searchIt.GetPixel(nPos) < (*flatPtr).second.bounds_min )
          { // If this is a boundary pixel && has a lesser value than
          // the currently recorded value...
          (*flatPtr).second.bounds_min = searchIt.GetPixel(nPos);
          (*flatPtr).second.min_label_ptr = labelIt[nPos];
          }
        if ( searchIt.GetPixel(nCenter) == searchIt.GetPixel(nPos) )
          {
          if ( labelIt.GetPixel(nPos) != NULL_LABEL )
            {
            // Pick up any equivalencies we missed before.
            equivalentLabels->Add(labelIt.GetPixel(nCenter),
                                  labelIt.GetPixel(nPos));
            }
          // If the following is encountered, it means that there is a
          // logic flaw in the first pass of this algorithm where flat
          // regions are initially detected and linked.
#ifndef NDEBUG
          else itkDebugMacro("An unexpected but non-fatal error has occurred.");
#endif
          }
              
        }
      }
    }

  // Merge the flat regions that we identified as connected components.
  Self::MergeFlatRegions(flatRegions, equivalentLabels);

  // Relabel the image with the merged regions.
  Self::RelabelImage(output, region, equivalentLabels);
}

template <class TInputImage>
void Segmenter<TInputImage>
::GradientDescent(InputImageTypePointer img,
                  ImageRegionType region)
{
  typename OutputImageType::Pointer output = this->GetOutputImage();

  InputPixelType minVal;
  unsigned int i, nPos;
  typename InputImageType::OffsetType  moveIndex;
  unsigned long newLabel;
  std::stack< unsigned long * > updateStack;
  
  //
  // Set up our iterators.
  //
  typename ConstNeighborhoodIterator<InputImageType>::RadiusType rad;
  typename NeighborhoodIterator<OutputImageType>::RadiusType zeroRad;
  for (i = 0; i < ImageDimension; ++i)
    {
    rad[i] = 1;
    zeroRad[i] = 0;
    }
  ConstNeighborhoodIterator<InputImageType>
    valueIt(rad, img, region);
  NeighborhoodIterator<OutputImageType>
    labelIt(zeroRad, output, region);
  ImageRegionIterator<OutputImageType> it(output, region);

  //
  // Sweep through the image and trace all unlabeled
  // pixels to a labeled region
  //
  for (it = it.Begin(); ! it.IsAtEnd(); ++it)
    {
    if ( it.Get() == NULL_LABEL )
      {
      valueIt.SetLocation(it.GetIndex());
      labelIt.SetLocation(it.GetIndex());
      newLabel = NULL_LABEL;               // Follow the path of steep-
      while( newLabel == NULL_LABEL )      // est descent until a label
        {                                  // is found.
        updateStack.push(labelIt.GetCenterPointer());
        minVal = valueIt.GetPixel(m_Connectivity.index[0]);
        moveIndex = m_Connectivity.direction[0];
        for (unsigned int ii = 1; ii < m_Connectivity.size; ++ii)
          {
          nPos = m_Connectivity.index[ii];
          if ( valueIt.GetPixel(nPos) < minVal) 
            {
            minVal = valueIt.GetPixel(nPos);
            moveIndex = m_Connectivity.direction[ii];
            }
          }
        valueIt += moveIndex;
        labelIt += moveIndex;
        newLabel = labelIt.GetPixel(0);
        }
        
      while( ! updateStack.empty() ) // Update all the pixels we've traversed
        {
        *(updateStack.top()) = newLabel;
        updateStack.pop();
        }
      }
    }
  
}

template <class TInputImage>
void Segmenter<TInputImage>
::DescendFlatRegions(flat_region_table_t &flatRegionTable,
                     ImageRegionType imageRegion)
{

  typename OutputImageType::Pointer output = this->GetOutputImage();
  // Assumes all pixels are labeled in the image.  Steps through the flat
  // regions and equates each one with the label at its lowest boundary
  // point. Flat basins are preserved as their own regions. The output image is
  // relabeled to reflect these equivalencies. 
  EquivalencyTable::Pointer equivalentLabels = EquivalencyTable::New();

  for (typename flat_region_table_t::const_iterator region = flatRegionTable.begin();
       region != flatRegionTable.end(); ++region)
    {
    if ( ((*region).second.bounds_min < (*region).second.value)
         && (! (*region).second.is_on_boundary) )
      {
      equivalentLabels->Add((*region).first, *((*region).second.min_label_ptr));
      }
    }
  
  equivalentLabels->Flatten();
  Self::RelabelImage(output, imageRegion, equivalentLabels);
}


template <class TInputImage>
void Segmenter<TInputImage>
::UpdateSegmentTable(InputImageTypePointer input, ImageRegionType region)
{
  edge_table_hash_t edgeHash;
  edge_table_t tempEdgeTable;

  typename edge_table_hash_t::iterator edge_table_entry_ptr;
  typename edge_table_t::iterator      edge_ptr;

  unsigned int i, nPos;
  typename NeighborhoodIterator<OutputImageType>::RadiusType hoodRadius;
  typename SegmentTableType::segment_t *segment_ptr;
  typename SegmentTableType::segment_t temp_segment;
  unsigned long segment_label;

  InputPixelType lowest_edge;
  
  // Grab the data we need.
  typename OutputImageType::Pointer output    = this->GetOutputImage();
  typename SegmentTableType::Pointer segments = this->GetSegmentTable();
 
  // Set up some iterators.
  for (i = 0; i < ImageDimension; i++)
    {
    hoodRadius[i] = 1;
    }
  ConstNeighborhoodIterator<InputImageType> searchIt(hoodRadius,input,region);
  NeighborhoodIterator<OutputImageType> labelIt(hoodRadius, output, region);

  unsigned long hoodCenter = searchIt.Size() >> 1;

  for (searchIt.GoToBegin(), labelIt.GoToBegin(); ! searchIt.IsAtEnd();
       ++searchIt, ++labelIt)
    {
    segment_label = labelIt.GetPixel(hoodCenter);
 
    // Find the segment corresponding to this label
    // and update its minimum value if necessary.
    segment_ptr = segments->Lookup(segment_label);
    edge_table_entry_ptr = edgeHash.find(segment_label);
    if (segment_ptr == 0) // This segment not yet identified.
      {                   // So add it to the table.
      temp_segment.min = searchIt.GetPixel(hoodCenter);
      segments->Add(segment_label, temp_segment);
      typedef typename edge_table_hash_t::value_type ValueType;
      edgeHash.insert(ValueType(segment_label,
                                      tempEdgeTable) );
          
      edge_table_entry_ptr = edgeHash.find(segment_label);
      }
    else if (searchIt.GetPixel(hoodCenter) < segment_ptr->min)
      {
      segment_ptr->min = searchIt.GetPixel(hoodCenter);
      }

    // Look up each neighboring segment in this segment's edge table.
    // If an edge exists, compare (and reset) the minimum edge value.
    // Note that edges are located *between* two adjacent pixels and
    // the value is taken to be the maximum of the two adjacent pixel
    // values.
    for (i = 0; i < m_Connectivity.size; ++i)
      {
      nPos = m_Connectivity.index[i];
      if (labelIt.GetPixel(nPos) != segment_label
          && labelIt.GetPixel(nPos) != NULL_LABEL)
        {
        if (searchIt.GetPixel(nPos) < searchIt.GetPixel(hoodCenter))
          lowest_edge = searchIt.GetPixel(hoodCenter); // We want the
        else lowest_edge = searchIt.GetPixel(nPos);    // max of the
        // adjacent pixels
              
        edge_ptr = (*edge_table_entry_ptr).second.find(labelIt.GetPixel(nPos));
        if ( edge_ptr == (*edge_table_entry_ptr).second.end() )
          {     // This edge has not been identified yet.
          typedef typename edge_table_t::value_type ValueType;
          (*edge_table_entry_ptr).second.insert(
            ValueType(labelIt.GetPixel(nPos), lowest_edge) ); 
          }
        else if (lowest_edge < (*edge_ptr).second)
          {
          (*edge_ptr).second = lowest_edge;
          }
        }
      }

    }

  //
  // Copy all of the edge tables into the edge lists of the
  // segment table.
  //
  unsigned long listsz;
  typename SegmentTableType::edge_list_t::iterator list_ptr;
  for (edge_table_entry_ptr = edgeHash.begin();
       edge_table_entry_ptr != edgeHash.end();
       edge_table_entry_ptr++)
    {
    // Lookup the corresponding segment entry
    segment_ptr = segments->Lookup((*edge_table_entry_ptr).first);
    if ( segment_ptr == 0 )
      {
      itkGenericExceptionMacro ( << "UpdateSegmentTable:: An unexpected and fatal error has occurred.");
      }

    // Copy into the segment list
    listsz = static_cast<unsigned long>( (*edge_table_entry_ptr).second.size() );
    segment_ptr->edge_list.resize(listsz);
    edge_ptr = (*edge_table_entry_ptr).second.begin();
    list_ptr = segment_ptr->edge_list.begin();
    while ( edge_ptr != (*edge_table_entry_ptr).second.end() )
      {
      list_ptr->label = (*edge_ptr).first;
      list_ptr->height= (*edge_ptr).second;
      edge_ptr++;
      list_ptr++;
      }

    // Clean up memory as we go
    (*edge_table_entry_ptr).second.clear();
    }
}

template <class TInputImage>
void Segmenter<TInputImage>
::BuildRetainingWall(InputImageTypePointer img,
                     ImageRegionType region,
                     InputPixelType value)
{
  unsigned int i;
  typename ImageRegionType::SizeType sz;
  typename ImageRegionType::IndexType idx;
  ImageRegionType reg;

  // Loop through the dimensions and populate the LOW and HIGH faces regions.
  for (i = 0; i < ImageDimension; ++i)
    {
    idx = region.GetIndex();       // LOW face
    sz  = region.GetSize();
    sz[i] = 1;
    reg.SetIndex(idx);
    reg.SetSize(sz);
    Segmenter::SetInputImageValues(img, reg, value);
    idx[i] = region.GetSize()[i] + region.GetIndex()[i] - 1;  // HIGH face
    reg.SetIndex(idx);
    Segmenter::SetInputImageValues(img, reg, value);
    }
}

/*
  ----------------------------------------------------------------------------
  Algorithm helper methods and debugging methods
  ----------------------------------------------------------------------------
*/
template <class TInputImage>
void Segmenter<TInputImage>
::SetInputImageValues(InputImageTypePointer img,
                 ImageRegionType region,
                 InputPixelType value)
{
  ImageRegionIterator<InputImageType> it(img, region);
  it = it.Begin();
  while (! it.IsAtEnd() )
    {
    it.Set(value);
    ++it;
    }
}


template <class TInputImage>
void Segmenter<TInputImage>
::SetOutputImageValues(OutputImageTypePointer img,
                 ImageRegionType region,
                 unsigned long value)
{
  ImageRegionIterator<OutputImageType> it(img, region);
  it = it.Begin();
  while (! it.IsAtEnd() )
    {
    it.Set(value);
    ++it;
    }
}

template <class TInputImage>
void Segmenter<TInputImage>
::MinMax(InputImageTypePointer img, ImageRegionType region,
         InputPixelType &min, InputPixelType &max)
{
  ImageRegionIterator<InputImageType> it(img, region);
  it = it.Begin();
  min = it.Value();
  max = it.Value();
  while (! it.IsAtEnd())
    {
    if (it.Get() > max) max = it.Get();
    if (it.Get() < min) min = it.Get();
    ++it;
    }
}

template <class TInputImage>
void Segmenter<TInputImage>
::MergeFlatRegions(flat_region_table_t &regions,
                   EquivalencyTable::Pointer eqTable)
{
  // Note that the labels must have no interdependencies.  That is,
  // every key must map to a value that is not itself a key in the
  // table. This means that you must always merge label->first with
  // label->second (a to b). EquivalencyTable can be converted to this
  // format with its Flatten() method.
  eqTable->Flatten();
  
  typename flat_region_table_t::iterator a, b;
  for (EquivalencyTable::ConstIterator it = eqTable->Begin();
       it != eqTable->End(); ++it)
    {
    if ( ((a = regions.find((*it).first)) == regions.end())
         || ((b = regions.find((*it).second)) == regions.end()) )
      {
      itkGenericExceptionMacro ( << "MergeFlatRegions:: An unexpected and fatal error has occurred.");
      }
      
    if ((*a).second.bounds_min < (*b).second.bounds_min)
      {
      (*b).second.bounds_min = (*a).second.bounds_min;
      (*b).second.min_label_ptr = (*a).second.min_label_ptr;
      }

    regions.erase(a);
    }
}

template <class TInputImage>
void Segmenter<TInputImage>
::RelabelImage(OutputImageTypePointer img,
               ImageRegionType region, 
               EquivalencyTable::Pointer eqTable)
{
  eqTable->Flatten();
  
  unsigned long temp;
  ImageRegionIterator<OutputImageType> it(img, region);

  it = it.Begin();
  while ( !it.IsAtEnd() )
    {
    temp = eqTable->Lookup(it.Get());
    if (temp != it.Get())  { it.Set(temp);}
    ++it; 
    }
}

template <class TInputImage>
void Segmenter<TInputImage>::
Threshold(InputImageTypePointer destination,
          InputImageTypePointer source,
          const ImageRegionType source_region,
          const ImageRegionType destination_region,
          InputPixelType threshold)
{
  ImageRegionIterator<InputImageType> dIt(destination, destination_region);
  ImageRegionIterator<InputImageType> sIt(source, source_region);

  dIt = dIt.Begin();
  sIt = sIt.Begin();

  // Assumes that source_region and destination region are the same size.  Does
  // no checking!!
  if (NumericTraits<InputPixelType>::is_integer)
    {
    // integral data type, if any pixel is at the maximum possible
    // value for the data type, then drop the value by one intensity
    // value. This the watershed algorithm to construct a "barrier" or
    // "wall" around the image that will stop the watershed without
    // requiring a expensive boundary condition checks.
    while ( ! dIt.IsAtEnd() )
      {
      InputPixelType tmp = sIt.Get();
      if ( tmp < threshold )
        {
        dIt.Set(threshold);
        }
      else if (tmp == NumericTraits<InputPixelType>::max())
        {
        dIt.Set(tmp - NumericTraits<InputPixelType>::One);
        }
      else
        {
        dIt.Set(tmp);
        }
      ++dIt;
      ++sIt;
      }
    }
  else
    {
    // floating point data, no need to worry about overflow
    while ( ! dIt.IsAtEnd() )
      {
      if ( sIt.Get() < threshold )
        {
        dIt.Set(threshold);
        }
      else
        {
        dIt.Set(sIt.Get());
        }
      ++dIt;
      ++sIt;
      }
    }
}

/*
  ----------------------------------------------------------------------------
  Pipeline methods
  ----------------------------------------------------------------------------
*/
template<class TInputImage>
typename Segmenter<TInputImage>::DataObjectPointer
Segmenter<TInputImage>
::MakeOutput(unsigned int idx)
{
  if (idx == 0)
    {  return static_cast<DataObject*>(OutputImageType::New().GetPointer());}
  else if (idx == 1)
    { return static_cast<DataObject*>(SegmentTableType::New().GetPointer());}
  else if (idx == 2)
    { return static_cast<DataObject*>(BoundaryType::New().GetPointer());}
  else return 0;
}

  
template <class TInputImage>
void 
Segmenter<TInputImage>::UpdateOutputInformation()
{
  unsigned int i;
  // call the superclass' implementation of this method
  Superclass::UpdateOutputInformation();
  
  // get pointers to the input and output
  typename InputImageType::Pointer  inputPtr = this->GetInputImage();
  typename OutputImageType::Pointer outputPtr = this->GetOutputImage();

  if ( !inputPtr || !outputPtr )
    {
    return;
    }
  // we need to compute the output spacing, the output image size, and the
  // output image start index
  const typename InputImageType::SizeType&   inputSize
    = inputPtr->GetLargestPossibleRegion().GetSize();
  const typename InputImageType::IndexType&  inputStartIndex
    = inputPtr->GetLargestPossibleRegion().GetIndex();
  
  typename OutputImageType::SizeType     outputSize;
  typename OutputImageType::IndexType    outputStartIndex;
  
  for (i = 0; i < OutputImageType::ImageDimension; i++)
    {
    outputSize[i] = inputSize[i];
    outputStartIndex[i]  = inputStartIndex[i];
    }

  typename OutputImageType::RegionType outputLargestPossibleRegion;
  outputLargestPossibleRegion.SetSize( outputSize );
  outputLargestPossibleRegion.SetIndex( outputStartIndex );
  
  outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
}

template <class TInputImage>
void Segmenter<TInputImage>::GenerateInputRequestedRegion()
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();
  
  // get pointers to the input and output
  typename InputImageType::Pointer  inputPtr  = this->GetInputImage();
  typename OutputImageType::Pointer outputPtr = this->GetOutputImage();
  
  if ( !inputPtr || !outputPtr )
    {
    return;
    }
  
  // 
  // FOR NOW WE'LL JUST SET THE INPUT REGION TO THE OUTPUT REGION
  // AND OVERRIDE THIS LATER
  //
  inputPtr->SetRequestedRegion( outputPtr->GetRequestedRegion() );
  
}

template <class TInputImage>
void 
Segmenter<TInputImage>
::GenerateOutputRequestedRegion(DataObject *output)
{
  // Only the Image output need to be propagated through.
  // No choice but to use RTTI here.
  ImageBase<ImageDimension> *imgData;
  ImageBase<ImageDimension> *op;
  imgData = dynamic_cast<ImageBase<ImageDimension> * >(output);
  typename TInputImage::RegionType c_reg;

  if (imgData)
    {
    std::vector<ProcessObject::DataObjectPointer>::size_type idx;
    for (idx = 0; idx < this->GetOutputs().size(); ++idx)
      {
      if (this->GetOutputs()[idx] && this->GetOutputs()[idx] != output)
        {
        op = dynamic_cast<ImageBase<ImageDimension>
          *>(this->GetOutputs()[idx].GetPointer());

        if (op) this->GetOutputs()[idx]->SetRequestedRegion(output);
        }
      }
    }
} 

template <class TInputImage>
Segmenter<TInputImage>
::Segmenter()
{
  m_Threshold = 0.0;
  m_MaximumFloodLevel = 1.0;
  m_CurrentLabel = 1;
  m_DoBoundaryAnalysis = false;
  m_SortEdgeLists = true;
  m_Connectivity.direction = 0;
  m_Connectivity.index = 0;
  typename OutputImageType::Pointer img
    = static_cast<OutputImageType*>(this->MakeOutput(0).GetPointer());
  typename SegmentTableType::Pointer st
    = static_cast<SegmentTableType*>(this->MakeOutput(1).GetPointer());
  typename BoundaryType::Pointer bd
    = static_cast<BoundaryType*>(this->MakeOutput(2).GetPointer());
  this->SetNumberOfRequiredOutputs(3);
  this->ProcessObject::SetNthOutput(0, img.GetPointer());
  this->ProcessObject::SetNthOutput(1, st.GetPointer());
  this->ProcessObject::SetNthOutput(2, bd.GetPointer());

  // Allocate memory for connectivity 
  m_Connectivity.size = 2 * ImageDimension;
  m_Connectivity.index = new unsigned int[m_Connectivity.size];
  m_Connectivity.direction
    = new typename InputImageType::OffsetType[m_Connectivity.size];
}

template<class TInputImage>
void 
Segmenter<TInputImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
  os << indent << "SortEdgeLists: " << m_SortEdgeLists << std::endl;
  os << indent << "DoBoundaryAnalysis: " << m_DoBoundaryAnalysis << std::endl;
  os << indent << "Threshold: " << m_Threshold << std::endl;
  os << indent << "MaximumFloodLevel: " << m_MaximumFloodLevel << std::endl;
  os << indent << "CurrentLabel: " << m_CurrentLabel << std::endl;  
}

}// end namespace watershed
}// end namespace itk

#endif