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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkKLMRegionGrowImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2004/10/23 12:48:34 $
  Version:   $Revision: 1.57 $

  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 _itkKLMRegionGrowImageFilter_txx
#define _itkKLMRegionGrowImageFilter_txx
#include "itkKLMRegionGrowImageFilter.h"

namespace itk
{


template<class TInputImage, class TOutputImage>
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::KLMRegionGrowImageFilter(void):
  m_MaximumLambda(1000),
  m_NumberOfRegions(0),
  m_InternalLambda(0),
  m_InitialNumberOfRegions(0),
  m_BorderCandidate(NULL),
  m_InitialRegionArea(0)
{
  m_InitialRegionMean.set_size( InputImageVectorDimension );
  m_InitialRegionMean.fill(0);
  this->SetMaximumNumberOfRegions( 2 );
}


template<class TInputImage, class TOutputImage>
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::~KLMRegionGrowImageFilter()
{
}


/*
 * PrintSelf
 */
template <class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);

  os << indent << "KLM Region grow segmentation object" << std::endl;
  os << indent << "KLM Region grow image filter object" << std::endl;
  os << indent << "Maximum value of lambda parameter: " << m_MaximumLambda << std::endl;
  os << indent << "Current internal value of lambda parameter: " << m_InternalLambda << std::endl;
  os << indent << "Initial number of regions: " << m_InitialNumberOfRegions << std::endl;
  os << indent << "Current number of regions: " << m_NumberOfRegions << std::endl;

} // end PrintSelf


/*
 * GenerateInputRequestedRegion method.
 */
template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::GenerateInputRequestedRegion()
{
  // This filter requires all of the input image to be in the buffer
  InputImagePointer inputPtr =
    const_cast< InputImageType * >( this->GetInput() );
  inputPtr->SetRequestedRegionToLargestPossibleRegion();
}


/*
 * EnlargeOutputRequestedRegion method.
 */
template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::EnlargeOutputRequestedRegion(
  DataObject *output )
{
  // This filter requires all of the output image to be in the buffer
  TOutputImage *imgData;
  imgData = dynamic_cast< TOutputImage* >( output );
  imgData->SetRequestedRegionToLargestPossibleRegion();
}


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::GenerateData()
{
  // Run the KLM algorithm
  this->ApplyRegionGrowImageFilter();

  // Set the output labelled and allocate the memory
  OutputImagePointer outputPtr = this->GetOutput();

  // Allocate the output buffer memory
  outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
  outputPtr->Allocate();

  GenerateOutputImage();

} // end GenerateData


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::GenerateOutputImage()
{
  InputImageConstPointer inputImage     = this->GetInput();
  InputImageSizeType     inputImageSize = inputImage->GetBufferedRegion().GetSize();

  GridSizeType gridSize = this->GetGridSize();

  InputImageSizeType numRegionsAlongDim;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    numRegionsAlongDim[idim] = inputImageSize[idim] / gridSize[idim];

  // Walk through each atomic block and get the approximation image.
  // This is needed as the region labels associated each atomic
  // block (as it was for during the initialization stage), have changed.
  // The new region label points to the region that has the up-to-date
  // region intensity. After the updated region intensity is found,
  // each pixel is updated with the mean approximation value.

  OutputImagePointer outputImage = this->GetOutput();
  typedef typename TOutputImage::RegionType OutputRegionType;
  OutputRegionType region;
  region.SetSize( gridSize ); // Constant grid size

  OutputImageIndexType     tmpIndex;
  tmpIndex.Fill(0);

  for ( unsigned int iregion = 0; iregion < m_InitialNumberOfRegions; iregion++ )
    {

    region.SetIndex( tmpIndex * gridSize );

    // Convert the mean region value to the correct output format

    MeanRegionIntensityType tmpMeanValue;
    OutputImageVectorType   outMeanValue;
    typedef typename OutputImagePixelType::ValueType OutputValueType;

    tmpMeanValue = m_RegionsPointer[iregion]->GetMeanRegionIntensity();

    for ( unsigned int ivecdim = 0; ivecdim < InputImageVectorDimension; ivecdim++ )
      outMeanValue[ivecdim] =
        static_cast< OutputValueType >( tmpMeanValue[ivecdim] );

    // Fill the region with the mean value

    OutputImageIterator outputIt( outputImage, region );
    OutputImageIterator outputItEnd = outputIt.End();

    while ( outputIt < outputItEnd )
      {
      outputIt.Set( outMeanValue );
      ++outputIt;
      }

    // Calculate next grid index
    long tmpVal = 1;
    for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
      {
      tmpIndex[ idim ]++;
      tmpVal *= numRegionsAlongDim[idim];
      if ( ( iregion + 1 ) % tmpVal != 0 ) break;
      tmpIndex[ idim ] = 0;
      }

    } // end iregion loop

} // end GenerateOutputImage()


template<class TInputImage, class TOutputImage>
typename KLMRegionGrowImageFilter<TInputImage,TOutputImage>::LabelImagePointer
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::GetLabelledImage()
{

  // Allocate the memory for the labelled image

  LabelImagePointer labelImagePtr = LabelImageType::New();

  typename LabelImageType::SizeType labelImageSize =
    this->GetInput()->GetBufferedRegion().GetSize();

  LabelImageIndexType labelImageIndex;
  labelImageIndex.Fill(0);

  typename LabelImageType::RegionType labelImageRegion;

  labelImageRegion.SetSize( labelImageSize );
  labelImageRegion.SetIndex( labelImageIndex );

  labelImagePtr->SetLargestPossibleRegion( labelImageRegion );
  labelImagePtr->SetBufferedRegion( labelImageRegion );
  labelImagePtr->Allocate();

  labelImagePtr = GenerateLabelledImage( labelImagePtr );

  return labelImagePtr;

} // end GetLabelledImage()


template<class TInputImage, class TOutputImage>
typename KLMRegionGrowImageFilter<TInputImage,TOutputImage>::LabelImagePointer
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::GenerateLabelledImage( LabelImageType *labelImagePtr )
{

  InputImageConstPointer inputImage     = this->GetInput();
  InputImageSizeType     inputImageSize = inputImage->GetBufferedRegion().GetSize();

  GridSizeType gridSize = this->GetGridSize();

  InputImageSizeType     numRegionsAlongDim;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    numRegionsAlongDim[idim] = inputImageSize[idim] / gridSize[idim];

  // Walk through each atomic block and get the new unique labels
  // representing the final segmentation.
  // This is needed as the region labels associated each atomic
  // block (as it was during the initialization stage), have changed.

  InputRegionType region;
  region.SetSize( gridSize ); // Constant grid size

  InputImageIndexType tmpIndex;
  tmpIndex.Fill(0);

  for ( unsigned int iregion = 0; iregion < m_InitialNumberOfRegions; iregion++ )
    {

    region.SetIndex( tmpIndex * gridSize );

    // Fill the region with the label

    RegionLabelType newRegionLabel = m_RegionsPointer[iregion]->GetRegionLabel();

    LabelImageIterator labelIt( labelImagePtr, region );
    LabelImageIterator labelItEnd = labelIt.End();

    while ( labelIt < labelItEnd )
      {
      labelIt.Set( newRegionLabel );
      ++labelIt;
      }

    // Calculate next grid index
    long tmpVal = 1;
    for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
      {
      tmpIndex[ idim ]++;
      tmpVal *= numRegionsAlongDim[idim];
      if ( ( iregion + 1 ) % tmpVal != 0 ) break;
      tmpIndex[ idim ] = 0;
      }
    }

  // Return the reference to the labelled image
  return labelImagePtr;

} // end GenerateLabelledImage()


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::ApplyRegionGrowImageFilter()
{
  this->ApplyKLM();
} // end ApplyRegionGrowImageFilter()


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::ApplyKLM()
{

  // Region merging based on the minimum value of the current
  // lambdas associated with the borders. The border with the
  // smallest lambda value will be taken out for region merging.
  // This growing process is repeated until the number of current
  // different regions is not bigger than the desired and the
  // current minimum scale parameter is not less than the desired scale.

  this->InitializeKLM();

  while (( m_NumberOfRegions > this->GetMaximumNumberOfRegions() ) &&
         ( m_InternalLambda  < m_MaximumLambda))
    {
    this->MergeRegions();
    }

  this->ResolveRegions();

} // end ApplyKLM()


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::InitializeKLM()
{

  // Maximum number of regions requested must be greater than 0
  if( this->GetMaximumNumberOfRegions() <= 1 )
    {
    itkExceptionMacro( << "Number of requested regions must be 2 or more" );
    }

  // This implementation requires the image dimensions to be
  // multiples of the user specified grid sizes.

  InputImageConstPointer inputImage     = this->GetInput();
  InputImageSizeType     inputImageSize = inputImage->GetBufferedRegion().GetSize();
  GridSizeType           gridSize       = this->GetGridSize();
  typename TInputImage::SpacingType
                         spacing        = inputImage->GetSpacing();

  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {
    if( gridSize[idim] == 0 ||
        inputImageSize[idim] % gridSize[idim] != 0 )
      {
      itkExceptionMacro( << "Invalid grid size" );
      }
    }

  // Determine the regions first and intialize them

  InputImageSizeType numRegionsAlongDim;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    numRegionsAlongDim[idim] = inputImageSize[idim] / gridSize[idim];

  // Calculate the initial number of regions

  m_InitialNumberOfRegions = 1;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {
    m_InitialNumberOfRegions *= numRegionsAlongDim[idim];
    }

  if( m_InitialNumberOfRegions < this->GetMaximumNumberOfRegions() )
    {
    itkWarningMacro( << "Number of initial image regions is less than requested: reduce granularity of the grid");
    }

  // Set current number of regions

  m_NumberOfRegions = m_InitialNumberOfRegions;

  // Allocate and intialize memory to the regions in initial image block

  m_RegionsPointer.resize( m_InitialNumberOfRegions );
  for ( unsigned int iregion = 0; iregion < m_InitialNumberOfRegions; iregion++ )
    {
    m_RegionsPointer[iregion] = KLMSegmentationRegion::New();
    }

  // Label each region

  InputRegionType  region;
  region.SetSize( gridSize ); // Constant grid size

  InputImageIndexType     tmpIndex;
  tmpIndex.Fill(0);

  for ( RegionLabelType iregion = 0; iregion < m_InitialNumberOfRegions; iregion++ )
    {

    region.SetIndex( tmpIndex * gridSize );

    InitializeRegionParameters( region );
    m_RegionsPointer[iregion]->SetRegionParameters( m_InitialRegionMean,
      m_InitialRegionArea, iregion + 1 );

    // Calculate next grid index
    long tmpVal = 1;
    for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
      {
      tmpIndex[ idim ]++;
      tmpVal *= numRegionsAlongDim[idim];
      if ( ( iregion + 1 ) % tmpVal != 0 ) break;
      tmpIndex[ idim ] = 0;
      }

    } // end iregion loop

  // Determine the borders next and intialize them

  // Allocate and intialize memory to the borders

  unsigned int numberOfBorders = 0;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {
    unsigned int tmpVal = 1;
    for ( unsigned int jdim = 0; jdim < InputImageDimension; jdim++ )
      {
      tmpVal *= ( jdim == idim ?
        numRegionsAlongDim[jdim] - 1 :
        numRegionsAlongDim[jdim] );
      }
    numberOfBorders += tmpVal;
    }

  // Allow a single region to pass through; this memory would not be
  // used but the memory allocation and free routine will throw
  // exception otherwise.
  if( numberOfBorders == 0 )
    {
    itkExceptionMacro( << "Number of initial regions must be 2 or more: reduce granularity of the grid" );
    }

  m_BordersPointer.resize( numberOfBorders );
  for ( unsigned int k = 0; k < m_BordersPointer.size(); k++ )
    {
    m_BordersPointer[k] = KLMSegmentationBorder::New();
    }

  /* the following initialization of the borders ensures that
     each border is assigned region1 and region2 such that,
     the label of region1 is less than the label of region2/
     and, that when a new border is added to a region,
     PushBack can be used for region1 and PushFront can be used
     for region2.  This will ensure that the borders are
     sorted by increased labels for region1 then region2 */

  m_TotalBorderLength = 0;
  unsigned int borderCounter   = 0;

  // Along each dimension, visit every border between two regions

  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {

    // along the dimension there is one less border than there are regions
    InputImageSizeType numBordersAlongDim = numRegionsAlongDim;
    numBordersAlongDim[idim]--;

    // compute number of borders to be seen this dimension and area of
    // each border
    unsigned int numBorderThisDim = 1;
    double       borderLengthTmp  = 1;
    for ( unsigned int jdim = 0; jdim < InputImageDimension; jdim++ )
      {
      numBorderThisDim *= numBordersAlongDim[jdim];
      borderLengthTmp  *= ( jdim == idim ? 1 : gridSize[ jdim ] * spacing[ jdim ] );
      }

    // index to atomic region1 and atomic region2
    InputImageIndexType  indexRegion1;
    InputImageIndexType  indexRegion2;
    indexRegion1.Fill(0);
    indexRegion2.Fill(0);
    indexRegion2[idim]++;

    for ( unsigned int iborder = 0; iborder < numBorderThisDim; iborder++ )
      {

      if ( borderCounter >= numberOfBorders )
        {
        itkExceptionMacro( << "KLM initialization is incorrect" );
        } // end if

      // Load the border of interest
      KLMSegmentationBorderPtr pcurrentBorder = m_BordersPointer[borderCounter];

      // Set the length of the border
      pcurrentBorder->SetBorderLength( borderLengthTmp );

      // m_TotalBorderLength is used as a sanity check
      m_TotalBorderLength += borderLengthTmp;

      // Find the two neighbor regions

      unsigned int intRegion1Index = 0;
      unsigned int intRegion2Index = 0;
      long tmpVal = 1;
      for ( unsigned int jdim = 0; jdim < InputImageDimension; jdim++ )
        {
        intRegion1Index += indexRegion1[jdim] * tmpVal;
        intRegion2Index += indexRegion2[jdim] * tmpVal;
        tmpVal *= numRegionsAlongDim[jdim];
        } // end jdim loop

      KLMSegmentationRegionPtr pRegion1 = m_RegionsPointer[ intRegion1Index ];
      KLMSegmentationRegionPtr pRegion2 = m_RegionsPointer[ intRegion2Index ];

      // Attach the region border off lesser label value to region1
      // Attach the region border of the greater label value to region2
      pcurrentBorder->SetRegion1( pRegion1 );
      pcurrentBorder->SetRegion2( pRegion2 );

      // The current border is linked to the region1 and region2
      // Initialize the border in the region objects
      pRegion1->PushBackRegionBorder( pcurrentBorder );
      pRegion2->PushFrontRegionBorder( pcurrentBorder );

      // Compute the scale parameter lambda
      pcurrentBorder->EvaluateLambda();

      // Increment the border counter to go to the next border
      borderCounter++;

      // Calculate next indices to atomic region1 and atomic region2
      tmpVal = 1;
      for ( unsigned int jdim = 0; jdim < InputImageDimension; jdim++ )
        {
        indexRegion1[ jdim ]++;
        tmpVal *= numBordersAlongDim[jdim];
        if ( ( iborder + 1 ) % tmpVal != 0 ) break;
        indexRegion1[ jdim ] = 0;
        }
      indexRegion2 = indexRegion1;
      indexRegion2[idim]++;

      }

    }

  // For DEBUG purposes
  if ( this->GetDebug() )
    {
    PrintAlgorithmRegionStats();
    PrintAlgorithmBorderStats();
    }

  // Verification of the initialization process

  double actualBorderLength = 0;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {
    double tmpDblVal = 1;
    for ( unsigned int jdim = 0; jdim < InputImageDimension; jdim++ )
      {
      tmpDblVal *= ( jdim == idim ?
        numRegionsAlongDim[jdim] - 1 : inputImageSize[jdim] * spacing[jdim] );
      }
    actualBorderLength += tmpDblVal;
    }

  if ( m_TotalBorderLength != actualBorderLength )
    {
    itkExceptionMacro( << "KLM initialization is incorrect" );
    } // end if
  else
    {
    itkDebugMacro( << "Passed initialization");
    } // end else

  // Allocate memory to store the array of pointers that point to the
  // static border objects

  m_BordersDynamicPointer.resize( m_BordersPointer.size() );
  for ( unsigned int k = 0; k < m_BordersDynamicPointer.size(); k++ )
    m_BordersDynamicPointer[ k ].m_Pointer = m_BordersPointer[k];

  // For DEBUG purposes
  if ( this->GetDebug() )
    {
    for ( unsigned int k = 0; k < m_BordersDynamicPointer.size(); k++ )
      itkDebugMacro( << m_BordersDynamicPointer[ k ].m_Pointer );
    }

  std::stable_sort(m_BordersDynamicPointer.begin(),
                  (m_BordersDynamicPointer.end()),
                  std::greater <KLMDynamicBorderArray<BorderType> >());

  m_BorderCandidate = &(m_BordersDynamicPointer[ m_BordersDynamicPointer.size() - 1 ]);
  m_InternalLambda = m_BorderCandidate->m_Pointer->GetLambda();

  if ( m_InternalLambda < 0.0 )
    itkExceptionMacro( << "KLM initialization is incorrect" );
} // end InitializeKLM()


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::InitializeRegionParameters( InputRegionType region )
{

  // Get a pointer to the image
  InputImageConstPointer   inputImage = this->GetInput();

  // Set the iterators and the pixel type definition for the input image
  InputImageConstIterator inputIt( inputImage, region );
  InputImageConstIterator inputItEnd = inputIt.End();

  // Variable to store the input pixel vector value
  InputImageVectorType     inputPixelVec;

  // Calculate V[0] for the constant model facet for the Region Growing
  // algorithm

  m_InitialRegionMean.fill(0);

  while ( inputIt < inputItEnd )
    {
    inputPixelVec = inputIt.Value();

    for ( unsigned int ivecdim = 0; ivecdim < InputImageVectorDimension; ivecdim++ )
      m_InitialRegionMean[ivecdim] += inputPixelVec[ivecdim];

    ++inputIt;
    }

  // Calculate the area and the mean associated with the region
  GridSizeType gridSize = this->GetGridSize();
  typename TInputImage::SpacingType spacing  = inputImage->GetSpacing();
  m_InitialRegionArea = 1;
  for ( unsigned int idim = 0; idim < InputImageDimension; idim++ )
    {
    m_InitialRegionArea *= gridSize[idim] * spacing[idim];
    }
  m_InitialRegionMean /= m_InitialRegionArea;

} // end InitializeRegionParameters


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::MergeRegions()
{
  itkDebugMacro( << "--------------------");
  itkDebugMacro( << "   Merging Regions  ");

  // Subtract border length before removing it
  m_TotalBorderLength -= m_BorderCandidate->m_Pointer->GetBorderLength();
  if ( m_TotalBorderLength <= 0 ) itkExceptionMacro( << "KLM algorithm error" );

  // Two regions are associated with the candidate border
  KLMSegmentationRegion *pRegion1;
  KLMSegmentationRegion *pRegion2;

  pRegion1 = m_BorderCandidate->m_Pointer->GetRegion1();
  pRegion2 = m_BorderCandidate->m_Pointer->GetRegion2();

  // For consistency, always assign smaller label: this affects
  // GenerateOutputImage and GenerateLabelledImage
  if( pRegion1->GetRegionLabel() >= pRegion2->GetRegionLabel() )
    {
    itkExceptionMacro( << "Invalid region labelling" );
    }

  // Add the new region's parameter data to the old.
  pRegion1->CombineRegionParameters( pRegion2 );

  // Remove the common region border from region 1 and region 2
  pRegion1->DeleteRegionBorder( m_BorderCandidate->m_Pointer );
  pRegion2->DeleteRegionBorder( m_BorderCandidate->m_Pointer );

  // Assign new equivalence label to the old region and update
  // the region borders, this is needed for ResolveRegions()
  pRegion2->ResetRegionLabelAndUpdateBorders( pRegion1 );

  // Merge the borders and region borders of two regions
  pRegion1->SpliceRegionBorders( pRegion2 );

  // Do not need the old region borders anymore
  pRegion2->DeleteAllRegionBorders();

  // Recompute the lambda's for all the borders of region1
  pRegion1->UpdateRegionBorderLambda();

  // Remove the common region border from list of sorted borders.
  // The BorderCandidate is always the last element.
  // Set the iterator to very last value and then erase that location
  m_BordersDynamicPointer.erase( m_BordersDynamicPointer.end() - 1 );

  // Decrement for the one deleted border and a deleted region
  m_NumberOfRegions--;
  if ( m_BordersDynamicPointer.empty() ) itkExceptionMacro( << "KLM algorithm error" );

  // For DEBUG purposes
  if ( this->GetDebug() )
    {
    itkDebugMacro( << "First Region ");
    pRegion1->PrintRegionInfo();
    itkDebugMacro( << "Second Region ");
    pRegion2->PrintRegionInfo();
    }

  // If any duplicate borders are found during SpliceRegionBorders,
  // lambda is set to -1.0, and pRegion1 and pRegion2 are set NULL
  // so that after this sort, the duplicate border will be the last
  // entry in m_BordersDynamicPointer

  // Resort the border list based on the lambda values
  std::stable_sort(m_BordersDynamicPointer.begin(),
                  (m_BordersDynamicPointer.end()),
                  std::greater <KLMDynamicBorderArray<BorderType> >());

  // Assign new BorderCandidate (it is always the last element).
  // Set Pointer to BorderCandidate to the last element
  m_BorderCandidate = &(m_BordersDynamicPointer[ m_BordersDynamicPointer.size() - 1 ]);
  m_InternalLambda = m_BorderCandidate->m_Pointer->GetLambda();

  // Remove any duplicate borders found during SpliceRegionBorders:
  // lambda = -1.0,  pRegion1 and pRegion2 = NULL
  while ( m_BorderCandidate->m_Pointer->GetRegion1() == NULL ||
          m_BorderCandidate->m_Pointer->GetRegion2() == NULL )
    {

    m_BordersDynamicPointer.erase( m_BordersDynamicPointer.end() - 1 );

    // Decrement for the one deleted border
    if ( m_BordersDynamicPointer.empty() ) itkExceptionMacro( << "KLM algorithm error" );

    m_BorderCandidate = &(m_BordersDynamicPointer[ m_BordersDynamicPointer.size() - 1 ]);
    m_InternalLambda = m_BorderCandidate->m_Pointer->GetLambda();
    }

} // end MergeRegions


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::ResolveRegions()
{
  InputImageConstPointer   inputImage = this->GetInput();

  // Scan through the region labels to establish the correspondence
  // between the final region (and label) and the initial regions.

  // Set up the unique label container class
  typedef std::vector<RegionLabelType>   UnsignedIntVectorType;
  UnsignedIntVectorType                  uniqueLabelsVec;

  // Resolve region labels to contain only unique labels.
  // Go backward from largest to smallest region label
  std::vector< KLMSegmentationRegionPtr >::reverse_iterator regionsPointerIt =
    m_RegionsPointer.rbegin();
  std::vector< KLMSegmentationRegionPtr >::reverse_iterator regionsPointerItEnd =
    m_RegionsPointer.rend();
  RegionLabelType iregion = m_InitialNumberOfRegions;
  while ( regionsPointerIt != regionsPointerItEnd )
    {
    RegionLabelType origLabel = iregion;
    iregion--;
    RegionLabelType currLabel = m_RegionsPointer[ iregion ]->GetRegionLabel();

    // Unresolved chain
    if ( currLabel != origLabel )
      {
      // Resolve a chain of equivalences by first finding the end of the chain
      RegionLabelType uniqLabel = origLabel;
      RegionLabelType tmpLabel  = currLabel;
      while ( uniqLabel != tmpLabel )
        {
        // In memory, the regions go from 0 to label-1. Hence the -1 offset
        uniqLabel = tmpLabel;
        tmpLabel =  m_RegionsPointer[ uniqLabel - 1 ]->GetRegionLabel();
        } // end of chain (while loop)

      // Then re-walk the chain to change the label of each chain
      // member to be the last one just retrieved (uniqLabel)
      while ( currLabel != origLabel )
        {
        m_RegionsPointer[ origLabel - 1 ]->SetRegionLabel( uniqLabel );
        origLabel  = currLabel;
        currLabel = m_RegionsPointer[ origLabel - 1 ]->GetRegionLabel();
        } // end of while ( currLabel != origLabel )

      } // end of the if condition for detecting unresolved chain

    else // The original label is unique, record it
      {
      uniqueLabelsVec.push_back( origLabel );
      }

    regionsPointerIt++;
    } // end of all regions

  // Sort the unique labels
  std::sort(uniqueLabelsVec.begin(), uniqueLabelsVec.end() );

  // Remap sorted unique labels to consecutive values

  UnsignedIntVectorType remapLabelsVec( m_InitialNumberOfRegions, 0 );

  UnsignedIntVectorType::iterator uniqueLabelsVecIterator;
  uniqueLabelsVecIterator = uniqueLabelsVec.begin();

  RegionLabelType newLabelValue = 1;

  while ( uniqueLabelsVecIterator != uniqueLabelsVec.end() )
    {
    remapLabelsVec[ *uniqueLabelsVecIterator - 1 ] = newLabelValue;
    uniqueLabelsVecIterator++;
    newLabelValue++;
    }

  // Assign new consecutive labels
  for ( unsigned int iregion = 0; iregion < m_InitialNumberOfRegions; iregion++ )
    {
    RegionLabelType labelValue = m_RegionsPointer[ iregion ]->GetRegionLabel();

    RegionLabelType newLabelValue = remapLabelsVec[ labelValue - 1 ];
    double newAreaValue = m_RegionsPointer[ labelValue - 1 ]->GetRegionArea();
    MeanRegionIntensityType newMeanValue =
      m_RegionsPointer[ labelValue - 1 ]->GetMeanRegionIntensity();

    m_RegionsPointer[iregion]->SetRegionParameters( newMeanValue,
                                                    newAreaValue,
                                                    newLabelValue );

    } // end looping through the regions

} // end ResolveRegions()


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::PrintAlgorithmRegionStats()
{
  // Print the stats associated with all the regions
  for ( unsigned int k = 0; k < m_InitialNumberOfRegions; k++ )
    {
    int i = m_RegionsPointer[ k ]->GetRegionBorderSize();
    if ( i > 0 )
      {
      std::cout << "Stats for Region No: "
                << m_RegionsPointer[ k ]->GetRegionLabel()
                << std::endl;
      m_RegionsPointer[ k ]->PrintRegionInfo();
      }
    } // end region printloop

} // end PrintAlgorithmRegionStats


template<class TInputImage, class TOutputImage>
void
KLMRegionGrowImageFilter<TInputImage,TOutputImage>
::PrintAlgorithmBorderStats()
{
  // Print the stats associated with all the regions
  for ( unsigned int k = 0; k < m_BordersDynamicPointer.size(); k++)
    {
    std::cout << "Stats for Border No: " << ( k + 1 ) << std::endl;
    m_BordersDynamicPointer[ k ].m_Pointer->PrintBorderInfo();
    } // end region printloop

} // end PrintAlgorithmBorderStats


} // namespace itk

#endif
