// =========================================================================
// @author Leonardo Florez Valencia
// @email florez-l@javeriana.edu.co
// =========================================================================
#ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__
#define __fpa__Common__SliceBySliceRandomWalker__hxx__

#include <vector>
#include <itkExtractImageFilter.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkJoinSeriesImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
SliceBySliceRandomWalker( )
  : Superclass( )
{
  ivqITKInputConfigureMacro( InputLabels, TLabels );
  ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
~SliceBySliceRandomWalker( )
{
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
GenerateData( )
{
  typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
  typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
  typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;

  // Preare I/O objects
  const TImage* input = this->GetInput( );
  const TLabels* labels = this->GetInputLabels( );
  const TScalarImage* vesselness = this->GetInputVesselness( );
  this->AllocateOutputs( );

  // Some values
  typename TImage::RegionType region = input->GetRequestedRegion( );
  typename TImage::SizeType regionIndex = region.GetIndex( );
  typename TImage::SizeType sliceSize = region.GetSize( );
  int numSlices = sliceSize[ 2 ];
  sliceSize[ 2 ] = 0;

  // Composite image
  typename TScalarImage::Pointer composite = TScalarImage::New( );
  composite->CopyInformation( vesselness );
  composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
  composite->Allocate( );
  this->_Composite( composite, labels, vesselness );

  // Extract slices
  std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
  std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
  std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
  for( int i = 0; i < numSlices; ++i )
  {
    // Prepare extraction region
    typename TImage::IndexType sliceIndex = regionIndex;
    sliceIndex[ 2 ] += i;
    typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );

    // Extract regions
    typename _TSliceImage::Pointer input_slice;
    typename _TSliceLabels::Pointer labels_slice;
    typename _TSliceScalarImage::Pointer composite_slice;
    this->_Slice( input_slice, input, sliceRegion );
    this->_Slice( labels_slice, labels, sliceRegion );
    this->_Slice( composite_slice, composite, sliceRegion );

    // Update vectors with each image
    data3D[ i ] = input_slice;
    binaryTree[ i ] = labels_slice;
    fusion[ i ] = composite_slice;

  } // rof

  // Random walker slice-by-slice
  this->_GoDown( binaryTree, data3D, fusion, numSlices );
  this->_GoUp( binaryTree, data3D, fusion, numSlices );

  // Join results
  typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
  typename _TJoin::Pointer join = _TJoin::New( );
  for( int i = 0; i < numSlices; ++i )
    join->SetInput( i, binaryTree[ i ] );
  join->Update( );
  this->GetOutput( )->Graft( join->GetOutput( ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
_Composite(
  typename TScalarImage::Pointer& composite,
  const TLabels* labels,
  const TScalarImage* vesselness
  )
{
  // Min-Max vesselness values
  typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
  typename _TMinMax::Pointer minMax = _TMinMax::New( );
  minMax->SetImage( vesselness );
  minMax->Compute( );
  TScalar maxVess = minMax->GetMaximum( );
  typename TScalarImage::RegionType region = labels->GetRequestedRegion( );

  // Composite image
  typename TScalarImage::Pointer composite = TScalarImage::New( );
  composite->CopyInformation( vesselness );
  composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
  composite->Allocate( );
  itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
  itk::ImageRegionConstIterator< TLabels > l( labels, region );
  itk::ImageRegionIterator< TScalarImage > c( composite, region );
  v.GoToBegin( );
  l.GoToBegin( );
  c.GoToBegin( );
  for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
    c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
template< class _TSlicePtr, class _TInput >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
_Slice(
  _TSlicePtr& slice,
  const _TInput* input,
  typename _TInput::RegionType region
  )
{
  typedef typename _TSlicePTr::ObjectType _TSlice;
  typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;

  typename _TExtract::Pointer extract = _TExtract::New( );
  extract->SetInput( input );
  extract->SetExtractionRegion( region );
  extract->SetDirectionCollapseToIdentity( );
  extract->Update( );
  slice = extract->GetOutput( );
  slice->DisconnectPipeline( );
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
template< class _TBinaryTree, class _TData3D, class _TFusion >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
_GoDown(
  _TBinaryTree& binaryTree,
  const _TData3D& data3D,
  const _TFusion& fusion,
  int numSlices
  )
{
  typedef typename _TBinaryTree::value_type     _TSliceBinaryPtr;
  typedef typename _TData3D::value_type         _TSliceData3DPtr;
  typedef typename _TFusion::value_type         _TSliceFusionPtr;
  typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
  typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
  typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;

  int z = -1;
  while( z < numSlices - 2 )
  {
    z++;
    _TSliceBinaryPtr tmp = binaryTree[ z ];
    _TSliceBinaryPtr next = binaryTree[ z  + 1];
    _TSliceData3DPtr data = data3D[ z + 1 ];
    _TSliceFusionPtr vess = fusion[ z + 1 ];

    // Create seeds
    _TSliceBinaryPtr seeds = _TSliceBinary::New( );
    seeds->CopyInformation( tmp );
    seeds->SetRequestedRegion( tmp->GetRequestedRegion( ) );
    seeds->Allocate( );
    seeds->FillBuffer( 0 );

    /* TODO


       SliceType::RegionType      sliceRegion = tmp->GetLargestPossibleRegion();
       SliceType::SizeType                size = sliceRegion.GetSize();
       SliceType::SpacingType     spacing = tmp->GetSpacing();
       SliceType::PointType       origin = tmp->GetOrigin();
       SliceType::IndexType       start = sliceRegion.GetIndex();

       SliceType::Pointer seedsObject = SliceType::New();
       seedsObject->SetRegions(sliceRegion);
       seedsObject->SetOrigin(origin);
       seedsObject->SetSpacing(spacing);
       seedsObject->Allocate();
       seedsObject->FillBuffer(0);

       SliceType::Pointer seedsBackground = SliceType::New();
       seedsBackground->SetRegions(sliceRegion);
       seedsBackground->SetOrigin(origin);
       seedsBackground->SetSpacing(spacing);
       seedsBackground->Allocate();
       seedsBackground->FillBuffer(0);

       typedef itk::ImageRegionIteratorWithIndex<SliceType> SliceIteratorType;
       SliceIteratorType itSliceMori(tmp, tmp->GetLargestPossibleRegion());

       int len = 0;
       int numO = 0;
       int numB = 0;

       itSliceMori.GoToBegin();

       while (!itSliceMori.IsAtEnd())
       {
       if (itSliceMori.Get() != 0)
       len++;

       typename TImage::PixelType vessVal = vess->GetPixel(itSliceMori.GetIndex());
       typename TImage::PixelType imgVal = data->GetPixel(itSliceMori.GetIndex());
       typename TImage::PixelType bw = tmp->GetPixel(itSliceMori.GetIndex());

       if (bw != 0 && vessVal > this->m_VessTh * this->m_vessMax / 100)
       {
       seedsObject->SetPixel(itSliceMori.GetIndex(), 1);
       numO++;
       }

       itSliceMori++;
       }

       if (len == 0)
       {
       std::cout << "Slice " << z << ": - brak pikseli zakwalifikowanych do drzewa" << std::endl;
       continue;
       }

       SliceIteratorType itSliceVess(vess, vess->GetLargestPossibleRegion());
       itSliceVess.GoToBegin();

       while (!itSliceVess.IsAtEnd())
       {
       if (itSliceVess.Get() < this->m_VessTh * this->m_vessMax / 100)
       {
       seedsBackground->SetPixel(itSliceVess.GetIndex(), 1);
       numB++;
       }
       itSliceVess++;

       }

       std::cout << "slice: " << z << "seeds obj: " << numO << " seeds bkg " << numB << std::endl;
       // przygotowania do Random Walkera

       typename TImage::PixelType  *bufferImg = 0;
       typename TImage::PixelType  *bufferSeedsObj = 0;
       typename TImage::PixelType  *bufferSeedsBkg = 0;

       int dims[2], dims2[2], dims3[2];
       double spac[2], orig[2], spac2[2], orig2[2], spac3[2], orig3[2];

       image2buffer(bufferImg, dims, spac, orig, 2, data);
       image2buffer(bufferSeedsObj, dims2, spac2, orig2, 2, seedsObject);
       image2buffer(bufferSeedsBkg, dims3, spac3, orig3, 2, seedsBackground);

       // TO DO HERE
       int width = dims[0];
       int height = dims[1];

       PixelType *mask = new PixelType[width*height];
       double *probs = new double[2 * width*height];

       int i, j, no_seeds = 0;

       for (i = 0; i < width*height; i++)

       if (bufferSeedsObj[i] == 1 || bufferSeedsBkg[i] == 1)
       no_seeds++;

       int *seed_indexes = new int[no_seeds];
       int *seed_labels = new int[no_seeds];

       for (j = 0, i = 0; i < width*height; i++)
       {
       if (bufferSeedsObj[i] == 1)
       {
       seed_indexes[j] = i + 1;//=i+1 (i)
       seed_labels[j] = this->m_replaceValue;
       j++;
       }
       else if (bufferSeedsBkg[i] == 1)
       {
       seed_indexes[j] = i + 1;//=i+1 (i)
       seed_labels[j] = 0;
       j++;
       }
       }

       call_walker(mask, probs, bufferImg, height, width, seed_indexes, no_seeds, seed_labels, 2, this->m_beta);

       delete[] bufferSeedsObj;
       delete[] bufferSeedsBkg;
       delete[] seed_indexes;
       delete[] seed_labels;

       SliceType::Pointer img2 = buffer2image(mask, dims, spac, orig, 2);

       SliceType::PixelType *probs_img = new PixelType[width*height];
       for (i = 0; i < width*height; i++)
       probs_img[i] = floor(255 * probs[i]);

       SliceType::Pointer img3 = buffer2image(probs_img, dims, spac, orig, 2);

       typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
       MaxFilterType::Pointer maxFilter = MaxFilterType::New();
       maxFilter->SetInput(0, img2);
       maxFilter->SetInput(1, next);
       maxFilter->Update();

       SliceType::Pointer sth = maxFilter->GetOutput();
       ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());

       delete[] probs;
    */
  }
}

// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalarImage >
template< class _TBinaryTree, class _TData3D, class _TFusion >
void
_GoUp(
  _TBinaryTree& binaryTree,
  const _TData3D& data3D,
  const _TFusion& fusion,
  int numSlices
  )
{
}

#endif // __fpa__Common__SliceBySliceRandomWalker__hxx__
// eof - $RCSfile$