#ifndef __itkSignedMaurerParallelDistanceMapImageFilter_txx
#define __itkSignedMaurerParallelDistanceMapImageFilter_txx

#include "itkSignedMaurerParallelDistanceMapImageFilter.h"
#include "itkImageSliceConstIteratorWithIndex.h"
#include "itkImageSliceIteratorWithIndex.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkAddImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
#include "itkBinaryThresholdImageFilter.h"

#include "vnl/vnl_math.h"

#define indexOfArray(x,y,z) ((z*m_Size[1]*m_Size[0])+(y*m_Size[0])+x)

/** Function to synchronize by waiting for all threads. */
namespace itk
{
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::WaitForAll ()
{
  m_Barrier->Wait();
}

/** Allocate data, do conversion to values of 1 and 0, erode the image,
 *  sum the images, and do the distance transform. */
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::GenerateData()
{
  // Allocate image.
  this->GetOutput()->SetRegions(this->GetInput()->GetRequestedRegion());
  this->GetOutput()->Allocate();

  // Used for anisotropic image dimensions.
  m_SpacingSquared = this->GetOutput()->GetSpacing();
  m_SpacingSquared[0] = m_SpacingSquared[0] * m_SpacingSquared[0];
  m_SpacingSquared[1] = m_SpacingSquared[1] * m_SpacingSquared[1];
  m_SpacingSquared[2] = m_SpacingSquared[2] * m_SpacingSquared[2];   

  typedef BinaryThresholdImageFilter<InputImageType, InputImageType> 
    BinaryFilterType; 

  typename BinaryFilterType::Pointer binaryFilter = BinaryFilterType::New();
 
  if(m_BinaryConvert)
    {
    binaryFilter->SetLowerThreshold(m_ForegroundValue);
    binaryFilter->SetUpperThreshold(m_ForegroundValue);
    binaryFilter->SetInsideValue(1);
    binaryFilter->SetOutsideValue(0);
    binaryFilter->SetNumberOfThreads(this->GetNumberOfThreads());
    binaryFilter->SetInput(this->GetInput());
    }

  if(m_Signed)
    {
    typedef BinaryBallStructuringElement< 
      InputPixelType, 
      InputImageDimension  > StructuringElementType;  

    typedef BinaryErodeImageFilter< 
      InputImageType, 
      InputImageType, 
      StructuringElementType >     ErodeFilterType; 

    typename ErodeFilterType::Pointer binaryEroder = ErodeFilterType::New();

    // Erode the input image by 1.
    StructuringElementType structuringElement;
    structuringElement.SetRadius(1);
    structuringElement.CreateStructuringElement();
    binaryEroder->SetKernel(structuringElement);
    binaryEroder->SetErodeValue(m_ForegroundValue);
    binaryEroder->SetNumberOfThreads(this->GetNumberOfThreads());
    if(m_BinaryConvert)
      {
      binaryEroder->SetInput(binaryFilter->GetOutput());
      }
    else
      {
      binaryEroder->SetInput(this->GetInput());
      }

    typedef AddImageFilter<InputImageType,InputImageType,InputImageType> 
      AddFilterType;

    typename AddFilterType::Pointer sum = AddFilterType::New();

    // Add the eroded input image to the uneroded one.
    // This makes interior elements 2, edge elements 1, 
    // and exterior elements 0.
    sum->SetInput1(this->GetInput());
    sum->SetInput2(binaryEroder->GetOutput());
    sum->SetNumberOfThreads(this->GetNumberOfThreads());
    sum->Update();
  
    // Variable for this sum image used for signed calculations.
    m_BaseImage = sum->GetOutput();
    }
  else
    {
    // Variable for the input image used for unsigned calculations.
    if(m_BinaryConvert)
      {
      m_BaseImage = binaryFilter->GetOutput();
      }
    else
      {
      m_BaseImage = this->GetInput();
      }
    }

  // Class variable used particularly for the indexing macro.
  m_Size = m_BaseImage->GetLargestPossibleRegion().GetSize();


  // Allocate two arrays for storing closest feature points.
  m_ClosestExterior = 
    static_cast<ClosestPoint*>
      (calloc(m_Size[0]*m_Size[1]*m_Size[2],sizeof(struct ClosestPoint)));
  m_ClosestInterior = 
    static_cast<ClosestPoint*>
      (calloc(m_Size[0]*m_Size[1]*m_Size[2],sizeof(struct ClosestPoint)));

  if (m_ClosestExterior == NULL || m_ClosestInterior == NULL)
  {
    itkExceptionMacro( 
     << "Unable to allocate space for closest feature array." );
  }

  if(!m_Signed)
  {
    // Only one call is needed if interior distances are unnecessary
    ComputeDistanceTransform(m_ClosestExterior, 1, 0);
  }
  else
  {
    // Two calls are needed if interior distances are necessary
    ComputeDistanceTransform(m_ClosestExterior, 1, 0);

    ComputeDistanceTransform(m_ClosestInterior, 0, 1);
  }

  // Free matrices after completion of processing.
  free(m_ClosestExterior);
  free(m_ClosestInterior);
}

/** Split the requested region for the specified thread on the specified
 *  dimension. Based on SplitRequestedRegion from the ImageSource class. */
template <class TInputImage, class TOutputImage>
unsigned int
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::SplitRequestedRegionOnNPlane(unsigned int i, unsigned int num, 
                               OutputImageRegionType& splitRegion,
                               unsigned int dimension)
{
  // Get the output pointer.
  OutputImageType * outputPtr = this->GetOutput();
  const typename TOutputImage::SizeType& requestedRegionSize
    = outputPtr->GetRequestedRegion().GetSize();

  unsigned int splitAxis=dimension;
  typename TOutputImage::IndexType splitIndex;
  typename TOutputImage::SizeType splitSize;

  // Initialize the splitRegion to the output requested region.
  splitRegion = outputPtr->GetRequestedRegion();
  splitIndex = splitRegion.GetIndex();
  splitSize = splitRegion.GetSize();

  // We do not split on the Y dimension.
  if (requestedRegionSize[splitAxis] == 1)
    {
    itkDebugMacro("  Cannot Split");
    return 1;
    }

  // Determine the actual number of pieces that will be generated.
  typename TOutputImage::SizeType::SizeValueType range = 
    requestedRegionSize[splitAxis];
  unsigned int valuesPerThread = 
    (unsigned int)::ceil(range/(double)num);
  unsigned int maxThreadIdUsed = 
    (unsigned int)::ceil(range/(double)valuesPerThread) - 1;

  // Split the region.
  if (i < maxThreadIdUsed)
    {
    splitIndex[splitAxis] += i*valuesPerThread;
    splitSize[splitAxis] = valuesPerThread;
    }
  if (i == maxThreadIdUsed)
    {
    splitIndex[splitAxis] += i*valuesPerThread;
    // The last thread needs to process the rest of the dimension being split.
    splitSize[splitAxis] = splitSize[splitAxis] - i*valuesPerThread;
    }

  // Set the split region index and size variables.
  splitRegion.SetIndex( splitIndex );
  splitRegion.SetSize( splitSize );

  itkDebugMacro("  Split Piece: " << splitRegion );

  return maxThreadIdUsed + 1;
}

/** Start the multithreaded distance transform. */
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter< TInputImage, TOutputImage >
::ComputeDistanceTransform(ClosestPoint *closest, 
                           InputPixelType featureVal, bool writeOutput)
{
  // Set thread struct values from parameters.
  MaurerThreadStruct str;
  str.Filter = this;
  str.Closest = closest;
  str.Write = writeOutput;
  str.Feature = featureVal;

  // Set up multithreading.
  this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads());
  m_Barrier->Initialize(this->GetMultiThreader()->GetNumberOfThreads());
  this->GetMultiThreader()->
    SetSingleMethod(this->ComputeDistanceTransformThreaderCallback, &str);

  // Start multithreading.
  this->GetMultiThreader()->SingleMethodExecute();
}

/** Set up thread information and region split and procede with 
 *  the threaded distance transform. */
template<class TInputImage, class TOutputImage>
ITK_THREAD_RETURN_TYPE
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::ComputeDistanceTransformThreaderCallback( void * arg )
{
  MaurerThreadStruct *str;

  unsigned int totalZ;
  unsigned int totalX;
  unsigned int threadID;
  unsigned int threadCount;

  // Pull information from threading struct.
  threadID = ((MultiThreader::ThreadInfoStruct *)(arg))->ThreadID;
  threadCount = ((MultiThreader::ThreadInfoStruct *)(arg))->NumberOfThreads;

  str = (MaurerThreadStruct *)
          (((MultiThreader::ThreadInfoStruct *)(arg))->UserData);

  OutputImageRegionType splitRegionZPlane;
  totalZ = str->Filter->
    SplitRequestedRegionOnNPlane(threadID, threadCount, splitRegionZPlane, 2);

  OutputImageRegionType splitRegionXPlane;
  totalX = str->Filter->
    SplitRequestedRegionOnNPlane(threadID, threadCount, splitRegionXPlane, 0);

  // Call threaded function.
  str->Filter->
    ThreadedComputeDistanceTransform(splitRegionZPlane, splitRegionXPlane, 
                                     threadID, totalZ, totalX, str->Closest, 
                                     str->Feature, str->Write);

  return ITK_THREAD_RETURN_VALUE;
}

/** Find the closest features in X, then Y, then Z, and do the 
 *  final distance transform calculation. Synchronize between the 
 *  x/y stage and the z one. */
template< class TInputImage, class TOutputImage >
void
SignedMaurerParallelDistanceMapImageFilter< TInputImage, TOutputImage >
::ThreadedComputeDistanceTransform(const OutputImageRegionType&
                                   outputZPlanesRegionForThread,
                                   const OutputImageRegionType&
                                   outputXPlanesRegionForThread,
                                   unsigned int threadID, unsigned int totalZ,
                                   unsigned int totalX,
                                   ClosestPoint *closest, 
                                   InputPixelType featureVal,
                                   bool writeOutput)
{
  // Only do if there are enough batches of Z planes for this thread.
  if(threadID<totalZ)
    {
    // Compute closest features in X.
    distanceForXPlane(closest, featureVal,
                      outputZPlanesRegionForThread, threadID);

    // Compute closest features in Y.
    distanceForYPlane(closest, outputZPlanesRegionForThread, threadID);
    }


  // Barrier to avoid conflicting accesses in the next step
  this->WaitForAll();


  // Only do if there are enough batches of X planes for this thread
  if(threadID<totalX)
    {
    // Compute closest features in Z.
    distanceForZPlane(closest, outputXPlanesRegionForThread, 
                      threadID, writeOutput);
    }
}

/** Compute the closest features in X by looping through each X line 
 *  and updating each element's closest feature. */
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::distanceForXPlane(ClosestPoint *closest, InputPixelType featureVal, 
                    const OutputImageRegionType& threadOutputRegion, 
                    unsigned int threadID)
{
  int numOfFeatures;
  InputIndexType start = threadOutputRegion.GetIndex();
  InputSizeType size = threadOutputRegion.GetSize();

  int *xFeatures = static_cast<int*>(calloc(size[0],sizeof(int)));
    
  int x = start[0];
  int y = start[1];
  int z = start[2];

  typedef ImageSliceConstIteratorWithIndex<InputImageType> ConstIteratorType;
  ConstIteratorType iter( m_BaseImage, threadOutputRegion );
  iter.SetFirstDirection(0);
  iter.SetSecondDirection(1);

  // Iterate over the region.
  while( z<(size[2]+start[2]) )
    {
    while( !iter.IsAtEndOfSlice() )
      {
      numOfFeatures = 0;

      // Iterate over the X line.
      while( !iter.IsAtEndOfLine() )
        {
        ClosestPoint v = {-1,-1,-1};

        // Initialize the closest feature array.
        closest[indexOfArray(x,y,z)] = v; 
        InputPixelType meshVal = iter.Get();
       
        // Find features based on featureVal.
        if((m_Signed && meshVal >= featureVal && meshVal < (featureVal+2)) ||
           (!m_Signed && meshVal == featureVal))
          {
          xFeatures[numOfFeatures] = x;
          numOfFeatures++;
          }
        ++iter;
        ++x;
        }
      x=start[0];
      
      // Assign closest feature values.
      for(int f(1); f <= numOfFeatures-1;++f)
        {
        double middle_x = (xFeatures[f-1] + xFeatures[f])/2;
            while(x <= middle_x)
            {
            ClosestPoint v = {xFeatures[f-1],y,z};
            int index = indexOfArray(x,y,z);
            closest[index] = v;

            x++;
            }
          }

        // Assign the last closest feature value.
        if(numOfFeatures > 0)
          {
          while (x < (start[0]+size[0]))
            {
            ClosestPoint v = {xFeatures[numOfFeatures-1],y,z};
            closest[indexOfArray(x,y,z)] = v;

            x++;
          }
        }
      x=0;
      ++y;
      iter.NextLine();
      }
    y=0;
    ++z;
    iter.NextSlice();
    }

  free(xFeatures);

}

/** Compute the closest features in Y by looping through each Y line and 
 *  updating each element's closest feature. */
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::distanceForYPlane(ClosestPoint* closest, 
                    const OutputImageRegionType& threadOutputRegion, 
                    unsigned int threadID)
{
  ClosestPoint thisClosest;
  int thisX(0);
  int thisY(0);
  int prevX(0);
  int prevY(0);
  int numOfFeatures;
  double yMiddle(0);

  InputIndexType start = threadOutputRegion.GetIndex();
  InputSizeType size = threadOutputRegion.GetSize();

  int *xFeatures = static_cast<int*>(calloc(size[1],sizeof(int)));
  int *yFeatures = static_cast<int*>(calloc(size[1],sizeof(int)));

  double *left_y = 
    static_cast<double*>(calloc(size[1],sizeof(double)));

  int previousFeature(0);

  // Iterate over the region.
  for (int z=start[2]; z < (size[2]+start[2]); ++z)
    {
    for (int x=start[0]; x < (size[0]+start[0]); ++x)
      {
      numOfFeatures = 0;
 
      // Iterate over the Y line.
      for (int y=start[1]; y < (size[1]+start[1]); ++y)
        {
        thisClosest = closest[indexOfArray(x,y,z)];
        
        // If thisClosest is not empty.
        if(thisClosest.x != -1 && thisClosest.y != -1 && thisClosest.z != -1)
          {
          thisX = thisClosest.x;
          thisY = thisClosest.y;
          previousFeature = numOfFeatures-1;
          yMiddle = 0;
          while (previousFeature >= 0)
            {
            // Retain only the valid previous features.
            prevX = xFeatures[previousFeature];
            prevY = yFeatures[previousFeature];
            
            double a;
            double b;
            double c;
            if(m_UseImageSpacing)
              {
              a = m_SpacingSquared[0] * 
                  (prevX - thisX) * ((2.0*x) - (prevX + thisX));
                b = 2.0 * m_SpacingSquared[1] * (thisY - prevY);
              }
            else
              {
              a = (prevX - thisX) * ((2.0*x) - (prevX + thisX));
              b = 2.0 * (thisY - prevY);
              }
            c = (prevY + thisY)/2.0;

            yMiddle = (a/b)+c;

            if(yMiddle <= left_y[previousFeature])
              {
              //Get rid of the previous feature.
              numOfFeatures--;
              previousFeature--;
              }
            else
              {
              previousFeature = -1; //Stop Checking
              }
            }
          
          // Add to features.
          if(yMiddle < (size[1]+start[1]))
            {
            xFeatures[numOfFeatures] = thisX;
            yFeatures[numOfFeatures] = thisY;
            left_y[numOfFeatures] = yMiddle;
            numOfFeatures++;
            }
          }
        }
      int y=start[1];

     // Assign the closest feature values.
     for(int f(1); f <= numOfFeatures-1;++f) 
        {
        while(y <= left_y[f])
          {
          ClosestPoint v = {xFeatures[f-1],yFeatures[f-1],z};
          closest[indexOfArray(x,y,z)] = v;

          y++;
          }
        }

      // Assign the last closest feature value.
      if(numOfFeatures > 0) 
        {
        while(y < (size[1]+start[1]))
          {
          ClosestPoint v =
            {xFeatures[numOfFeatures-1],yFeatures[numOfFeatures-1],z};

          closest[indexOfArray(x,y,z)] = v;

          y++;
          }
        }
      }
    }

  free(xFeatures);
  free(yFeatures);
  free(left_y);

}

/** Compute the closest features in Z by looping through each Z line and 
 *  updating each element's closest feature.  To avoid repeating loops, do 
 *  final distance calculation here as well. */
template<class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::distanceForZPlane(ClosestPoint* closest, 
                    const OutputImageRegionType& threadOutputRegion, 
                    unsigned int threadID, bool writeOutput)
{
  ClosestPoint thisClosest;
  int thisX(0);
  int thisY(0);
  int thisZ(0);
  int numOfFeatures;
  double zMiddle(0);
  int previousFeature(0);

  typedef ImageSliceIteratorWithIndex<OutputImageType> IteratorType;
  IteratorType imageIter( this->GetOutput(), threadOutputRegion );
  imageIter.SetFirstDirection(2);
  imageIter.SetSecondDirection(1);

  typedef ImageSliceConstIteratorWithIndex<InputImageType> ConstIteratorType;
  ConstIteratorType meshIter( m_BaseImage, threadOutputRegion );
  meshIter.SetFirstDirection(2);
  meshIter.SetSecondDirection(1);

  InputIndexType start = threadOutputRegion.GetIndex();
  InputSizeType size = threadOutputRegion.GetSize();


  int *xFeatures = static_cast<int*>(calloc(size[2],sizeof(int)));
  int *yFeatures = static_cast<int*>(calloc(size[2],sizeof(int)));
  int *zFeatures = static_cast<int*>(calloc(size[2],sizeof(int)));
  double *zLeft = 
    static_cast<double*>(
      calloc(size[2],sizeof(double)));

  // Iterate over the region.
  for(int x=start[0]; x < (size[0]+start[0]); 
     ++x, imageIter.NextSlice(), meshIter.NextSlice())
    {
    for(int y=start[1]; y < (size[1]+start[1]); 
        ++y, imageIter.NextLine(), meshIter.NextLine()) 
      {
      numOfFeatures=0;

      // Iterate over the Z line.
      for (int z=start[2]; z < (size[2]+start[2]); ++z) 
        {
        // Get the closest feature from the Z plane.
        thisClosest = closest[indexOfArray(x,y,z)]; 

        // If thisClosest is not empty.
        if(thisClosest.x != -1 && thisClosest.y != -1 && thisClosest.z != -1)
          {
          thisX = thisClosest.x;
          thisY = thisClosest.y;
          thisZ = thisClosest.z;
          previousFeature = numOfFeatures-1;
          zMiddle = start[2];
          
          // Retain only valid features.
          while(previousFeature >= 0)
            {
            double a1;
            double a2;
            double b;
            double c;
            if(m_UseImageSpacing)
              {
              a1 =  m_SpacingSquared[0] * 
                    (xFeatures[previousFeature] - thisX) * 
                    ((2.0*x) - (xFeatures[previousFeature] + thisX));
              a2 = m_SpacingSquared[1] * (yFeatures[previousFeature] - thisY) * 
                    ((2.0*y) - (yFeatures[previousFeature] + thisY));
              b = 2.0 * m_SpacingSquared[2] * (thisZ - zFeatures[previousFeature]);
              }
            else
              {
              a1 = (xFeatures[previousFeature] - thisX) * 
                   ((2.0*x) - (xFeatures[previousFeature] + thisX));
              a2 = (yFeatures[previousFeature] - thisY) * 
                   ((2.0*y) - (yFeatures[previousFeature] + thisY));
              b = 2.0 * (thisZ - zFeatures[previousFeature]);
              }
            c = (zFeatures[previousFeature] + thisZ)/2.0;

            zMiddle = ((a1 + a2)/static_cast<double>(b)) + c;

            if(zMiddle <= zLeft[previousFeature])
              {
              numOfFeatures--;
              previousFeature--;
              }
            else
              {
              previousFeature = -1;
              }
            }

          // Add to features.
          if(zMiddle < (start[2]+size[2]))
            {
            xFeatures[numOfFeatures] = thisX;
            yFeatures[numOfFeatures] = thisY;
            zFeatures[numOfFeatures] = thisZ;
            zLeft[numOfFeatures] = zMiddle;
            numOfFeatures++;
            }
          } 
        
        } 
      
        
      int z=start[2];
      // Assign closest features or output.
      for (int f(1); f <= numOfFeatures-1; ++f) 
        {
        while(z <= zLeft[f])
          {
          ClosestPoint v = {xFeatures[f-1],yFeatures[f-1],zFeatures[f-1]};
          closest[indexOfArray(x,y,z)]=v;

          if(writeOutput)
            {
            OutputPixelType square = 0;
            v=m_ClosestExterior[indexOfArray(x,y,z)];

            OutputPixelType distance;
            if(m_UseImageSpacing)
              {
              distance = static_cast<OutputPixelType>(
                m_SpacingSquared[0]*(x-v.x)*(x-v.x) +
                m_SpacingSquared[1]*(y-v.y)*(y-v.y) +
                m_SpacingSquared[2]*(z-v.z)*(z-v.z));
              }
            else
              {
              distance =
                (x-v.x)*(x-v.x) +
                (y-v.y)*(y-v.y) +
                (z-v.z)*(z-v.z);
              }

            if(m_Signed)
              {
              v=m_ClosestInterior[indexOfArray(x,y,z)];
              OutputPixelType signedDistance;
              if(m_UseImageSpacing)
                {
                signedDistance = static_cast<OutputPixelType>(
                  m_SpacingSquared[0]*(x-v.x)*(x-v.x) +
                  m_SpacingSquared[1]*(y-v.y)*(y-v.y) +
                  m_SpacingSquared[2]*(z-v.z)*(z-v.z));
                }
              else
                {
                signedDistance =
                  (x-v.x)*(x-v.x) +
                  (y-v.y)*(y-v.y) +
                  (z-v.z)*(z-v.z);
                }

              if(meshIter.Get() == 2)
                {
                if(m_SquaredDistance)
                  {
                  square = signedDistance;
                  }
                else
                  {
                  square = static_cast<OutputPixelType>(sqrt(signedDistance));
                  }
                if(!m_InsideIsPositive)
                  {
                  square*=-1;
                  }
                }
              else
                {
                if(m_SquaredDistance)
                  {
                  square = distance;
                  }
                else
                  {
                  square = static_cast<OutputPixelType>(sqrt(distance));
                  }
                if(m_InsideIsPositive)
                  {
                  square*=-1;
                  }
                }

              imageIter.Set(square);
              }
            else
              {
              if(m_SquaredDistance)
                {
                square = distance;
                }
              else
                {
                square = static_cast<OutputPixelType>(sqrt(distance));
                }

              imageIter.Set(square);
              }

            ++meshIter;
            ++imageIter;
            }
          z++; 
          }
        }

      // Assign last closest feature or output.
      if(numOfFeatures > 0)
        {
        while(z < (size[2]+start[2]))
          {
          ClosestPoint v = 
            {xFeatures[numOfFeatures-1],
             yFeatures[numOfFeatures-1],
             zFeatures[numOfFeatures-1]};
          closest[indexOfArray(x,y,z)]=v;
          if(writeOutput)
            {
            OutputPixelType square = 0;
            v=m_ClosestExterior[indexOfArray(x,y,z)];

            OutputPixelType distance;
            if(m_UseImageSpacing)
              {
              distance = static_cast<OutputPixelType>(
                m_SpacingSquared[0]*(x-v.x)*(x-v.x) +
                m_SpacingSquared[1]*(y-v.y)*(y-v.y) +
                m_SpacingSquared[2]*(z-v.z)*(z-v.z));
              }
            else
              {
              distance =
                (x-v.x)*(x-v.x) +
                (y-v.y)*(y-v.y) +
                (z-v.z)*(z-v.z);
              }

            if(m_Signed)
              {
              v=m_ClosestInterior[indexOfArray(x,y,z)];
                
              OutputPixelType signedDistance;
              if(m_UseImageSpacing)
                {
                signedDistance = static_cast<OutputPixelType>(
                  m_SpacingSquared[0]*(x-v.x)*(x-v.x) +
                  m_SpacingSquared[1]*(y-v.y)*(y-v.y) +
                  m_SpacingSquared[2]*(z-v.z)*(z-v.z));
                }
              else
                {
                signedDistance =
                  (x-v.x)*(x-v.x) +
                  (y-v.y)*(y-v.y) +
                  (z-v.z)*(z-v.z);
                }
                
              if(meshIter.Get() == 2)
                {

                if(m_SquaredDistance)
                  {
                  square = signedDistance;
                  }
                else
                  {
                  square = static_cast<OutputPixelType>(sqrt(signedDistance));
                  }
                if(!m_InsideIsPositive)
                  {
                  square*=-1;
                  }
                }
              else
                {
                if(m_SquaredDistance)
                  {
                  square = distance;
                  }
                else
                  {
                  square = static_cast<OutputPixelType>(sqrt(distance));
                  }

                if(m_InsideIsPositive)
                  {
                  square*=-1;
                  }
                }

              imageIter.Set(square);
              }
            else
              {
              if(m_SquaredDistance)
                {
                square = distance;
                }
              else
                {
                square = static_cast<OutputPixelType>(sqrt(distance));
                }

              imageIter.Set(square);
              }

            ++meshIter;
            ++imageIter;
            }

          z++;
          }
        }
      }
    }

  free(xFeatures);
  free(yFeatures);
  free(zFeatures);
  free(zLeft);
}

/**
 * Standard PrintSelf method
 */
template <class TInputImage, class TOutputImage>
void
SignedMaurerParallelDistanceMapImageFilter<TInputImage, TOutputImage>
::PrintSelf(
  std::ostream& os, 
  Indent indent) const
{
  Superclass::PrintSelf( os, indent );
  os << indent << "Foreground Value: " << m_ForegroundValue << std::endl;
  os << indent << "Spacing Squared: " << m_SpacingSquared << std::endl;
  os << indent << "Size: " << m_Size << std::endl;
  os << indent << "Inside Positive: " << m_InsideIsPositive << std::endl;
  os << indent << "Squared Distance: " << m_SquaredDistance << std::endl;
  os << indent << "Use Image Spacing: " << m_UseImageSpacing << std::endl;
  os << indent << "Convert to Binary Form: " << m_BinaryConvert << std::endl;
  os << indent << "Signed: " << m_Signed << std::endl;
  os << indent << "Barrier: " << m_Barrier << std::endl;
  os << indent << "Base Image: " << m_BaseImage << std::endl;
  os << indent << "Closest Interior Vortexes: " << m_ClosestInterior 
     << std::endl;
  os << indent << "Closest Exterior Vortexes: " << m_ClosestExterior 
     << std::endl;
}

} // End namespace itk.

#endif