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

  Program:   Visualization Toolkit
  Module:    vtkImageThresholdConnectivity.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/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 notice for more information.

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

#include "vtkImageThresholdConnectivity.h"

#include "vtkImageData.h"
#include "vtkImageIterator.h"
#include "vtkImageStencilData.h"
#include "vtkImageStencilIterator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTemplateAliasMacro.h"

#include <stack>

vtkStandardNewMacro(vtkImageThresholdConnectivity);
vtkCxxSetObjectMacro(vtkImageThresholdConnectivity, SeedPoints, vtkPoints);

//------------------------------------------------------------------------------
// Constructor sets default values
vtkImageThresholdConnectivity::vtkImageThresholdConnectivity()
{
  this->UpperThreshold = VTK_FLOAT_MAX;
  this->LowerThreshold = -VTK_FLOAT_MAX;
  this->SeedPoints = nullptr;
  this->ReplaceIn = 0;
  this->InValue = 0.0;
  this->ReplaceOut = 0;
  this->OutValue = 0.0;

  this->NeighborhoodRadius[0] = 0.0;
  this->NeighborhoodRadius[1] = 0.0;
  this->NeighborhoodRadius[2] = 0.0;
  this->NeighborhoodFraction = 0.5;

  this->SliceRangeX[0] = -VTK_INT_MAX;
  this->SliceRangeX[1] = VTK_INT_MAX;
  this->SliceRangeY[0] = -VTK_INT_MAX;
  this->SliceRangeY[1] = VTK_INT_MAX;
  this->SliceRangeZ[0] = -VTK_INT_MAX;
  this->SliceRangeZ[1] = VTK_INT_MAX;

  this->ActiveComponent = -1;

  this->ImageMask = vtkImageData::New();

  this->NumberOfInVoxels = 0;

  this->SetNumberOfInputPorts(2);
}

//------------------------------------------------------------------------------
vtkImageThresholdConnectivity::~vtkImageThresholdConnectivity()
{
  if (this->SeedPoints)
  {
    this->SeedPoints->Delete();
  }
  this->ImageMask->Delete();
}

//------------------------------------------------------------------------------
void vtkImageThresholdConnectivity::SetInValue(double val)
{
  if (val != this->InValue || this->ReplaceIn != 1)
  {
    this->InValue = val;
    this->ReplaceIn = 1;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
void vtkImageThresholdConnectivity::SetOutValue(double val)
{
  if (val != this->OutValue || this->ReplaceOut != 1)
  {
    this->OutValue = val;
    this->ReplaceOut = 1;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
// The values greater than or equal to the value match.
void vtkImageThresholdConnectivity::ThresholdByUpper(double thresh)
{
  if (this->LowerThreshold != thresh || this->UpperThreshold < VTK_FLOAT_MAX)
  {
    this->LowerThreshold = thresh;
    this->UpperThreshold = VTK_FLOAT_MAX;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
// The values less than or equal to the value match.
void vtkImageThresholdConnectivity::ThresholdByLower(double thresh)
{
  if (this->UpperThreshold != thresh || this->LowerThreshold > -VTK_FLOAT_MAX)
  {
    this->UpperThreshold = thresh;
    this->LowerThreshold = -VTK_FLOAT_MAX;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
// The values in a range (inclusive) match
void vtkImageThresholdConnectivity::ThresholdBetween(double lower, double upper)
{
  if (this->LowerThreshold != lower || this->UpperThreshold != upper)
  {
    this->LowerThreshold = lower;
    this->UpperThreshold = upper;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
int vtkImageThresholdConnectivity::FillInputPortInformation(int port, vtkInformation* info)
{
  if (port == 1)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageStencilData");
    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
  }
  else
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
  }
  return 1;
}

//------------------------------------------------------------------------------
void vtkImageThresholdConnectivity::SetStencilData(vtkImageStencilData* stencil)
{
  this->SetInputData(1, stencil);
}

//------------------------------------------------------------------------------
vtkImageStencilData* vtkImageThresholdConnectivity::GetStencil()
{
  if (this->GetNumberOfInputConnections(1) < 1)
  {
    return nullptr;
  }

  return vtkImageStencilData::SafeDownCast(this->GetExecutive()->GetInputData(1, 0));
}

//------------------------------------------------------------------------------
vtkMTimeType vtkImageThresholdConnectivity::GetMTime()
{
  vtkMTimeType mTime = this->MTime.GetMTime();
  vtkMTimeType pointsMTime;

  if (this->SeedPoints)
  {
    pointsMTime = this->SeedPoints->GetMTime();
    mTime = (pointsMTime > mTime ? pointsMTime : mTime);
  }

  return mTime;
}

//------------------------------------------------------------------------------
// seed struct: just a set of indices
class vtkFloodFillSeed
{
public:
  vtkFloodFillSeed()
  {
    store[0] = 0;
    store[1] = 0;
    store[2] = 0;
  }
  vtkFloodFillSeed(int i, int j, int k)
  {
    store[0] = i;
    store[1] = j;
    store[2] = k;
  }
  vtkFloodFillSeed(const vtkFloodFillSeed& seed)
  {
    store[0] = seed.store[0];
    store[1] = seed.store[1];
    store[2] = seed.store[2];
  }
  const int& operator[](int i) const { return store[i]; }
  vtkFloodFillSeed& operator=(const vtkFloodFillSeed& seed)
  {
    if (this == &seed)
    {
      return *this;
    }

    store[0] = seed.store[0];
    store[1] = seed.store[1];
    store[2] = seed.store[2];
    return *this;
  }

private:
  int store[3];
};

//------------------------------------------------------------------------------
// Make sure the thresholds are valid for the input scalar range
template <class IT>
void vtkImageThresholdConnectivityThresholds(
  vtkImageThresholdConnectivity* self, vtkImageData* inData, IT& lowerThreshold, IT& upperThreshold)
{
  if (self->GetLowerThreshold() < inData->GetScalarTypeMin())
  {
    lowerThreshold = static_cast<IT>(inData->GetScalarTypeMin());
  }
  else
  {
    if (self->GetLowerThreshold() > inData->GetScalarTypeMax())
    {
      lowerThreshold = static_cast<IT>(inData->GetScalarTypeMax());
    }
    else
    {
      lowerThreshold = static_cast<IT>(self->GetLowerThreshold());
    }
  }
  if (self->GetUpperThreshold() > inData->GetScalarTypeMax())
  {
    upperThreshold = static_cast<IT>(inData->GetScalarTypeMax());
  }
  else
  {
    if (self->GetUpperThreshold() < inData->GetScalarTypeMin())
    {
      upperThreshold = static_cast<IT>(inData->GetScalarTypeMin());
    }
    else
    {
      upperThreshold = static_cast<IT>(self->GetUpperThreshold());
    }
  }
}

//------------------------------------------------------------------------------
// Make sure the replacement values are within the output scalar range
template <class OT>
void vtkImageThresholdConnectivityValues(
  vtkImageThresholdConnectivity* self, vtkImageData* outData, OT& inValue, OT& outValue)
{
  if (self->GetInValue() < outData->GetScalarTypeMin())
  {
    inValue = static_cast<OT>(outData->GetScalarTypeMin());
  }
  else
  {
    if (self->GetInValue() > outData->GetScalarTypeMax())
    {
      inValue = static_cast<OT>(outData->GetScalarTypeMax());
    }
    else
    {
      inValue = static_cast<OT>(self->GetInValue());
    }
  }
  if (self->GetOutValue() > outData->GetScalarTypeMax())
  {
    outValue = static_cast<OT>(outData->GetScalarTypeMax());
  }
  else
  {
    if (self->GetOutValue() < outData->GetScalarTypeMin())
    {
      outValue = static_cast<OT>(outData->GetScalarTypeMin());
    }
    else
    {
      outValue = static_cast<OT>(self->GetOutValue());
    }
  }
}

//------------------------------------------------------------------------------
static void vtkImageThresholdConnectivityApplyStencil(
  vtkImageData* maskData, vtkImageStencilData* stencil, int extent[6])
{
  vtkImageStencilIterator<unsigned char> iter(maskData, stencil, extent);
  while (!iter.IsAtEnd())
  {
    unsigned char* beginptr = iter.BeginSpan();
    unsigned char* endptr = iter.EndSpan();
    unsigned char val = (iter.IsInStencil() ? 0 : 1);

    for (unsigned char* ptr = beginptr; ptr < endptr; ptr++)
    {
      *ptr = val;
    }

    iter.NextSpan();
  }
}

//------------------------------------------------------------------------------
// This templated function executes the filter for any type of data.
template <class IT, class OT>
void vtkImageThresholdConnectivityExecute(vtkImageThresholdConnectivity* self, vtkImageData* inData,
  vtkImageData* outData, vtkImageStencilData* stencil, vtkImageData* maskData, int outExt[6],
  int id, IT* inPtr, OT* outPtr, int& voxelCount)
{
  // Get active component (only one component is thresholded)
  int nComponents = outData->GetNumberOfScalarComponents();
  int activeComponent = self->GetActiveComponent();
  if (activeComponent < 0)
  {
    activeComponent = 0;
  }
  activeComponent = activeComponent % nComponents;

  // Get thresholds as input data type
  IT lowerThreshold, upperThreshold;
  vtkImageThresholdConnectivityThresholds(self, inData, lowerThreshold, upperThreshold);

  // Get replace values as output data type
  bool replaceIn = (self->GetReplaceIn() != 0);
  bool replaceOut = (self->GetReplaceOut() != 0);
  OT inValue, outValue;
  vtkImageThresholdConnectivityValues(self, outData, inValue, outValue);

  // Set the "outside" with either the input or the OutValue
  vtkImageIterator<IT> inIt(inData, outExt);
  vtkImageIterator<OT> outIt(outData, outExt);
  while (!outIt.IsAtEnd())
  {
    IT* inSI = inIt.BeginSpan();
    OT* outSI = outIt.BeginSpan();
    OT* outSIEnd = outIt.EndSpan();

    if (replaceOut)
    {
      if (nComponents == 1)
      {
        while (outSI < outSIEnd)
        {
          *outSI++ = outValue;
        }
      }
      else
      {
        // only color the active component, copy the rest
        while (outSI < outSIEnd)
        {
          int jj = 0;
          while (jj < activeComponent)
          {
            *outSI++ = static_cast<OT>(*inSI++);
            jj++;
          }
          *outSI++ = outValue;
          inSI++;
          jj++;
          while (jj < nComponents)
          {
            *outSI++ = static_cast<OT>(*inSI++);
            jj++;
          }
        }
      }
    }
    else
    {
      while (outSI < outSIEnd)
      {
        *outSI++ = static_cast<OT>(*inSI++);
      }
    }
    inIt.NextSpan();
    outIt.NextSpan();
  }

  // Get the extent for the flood fill, and clip with the input extent
  int extent[6];
  int inExt[6];
  self->GetSliceRangeX(&extent[0]);
  self->GetSliceRangeY(&extent[2]);
  self->GetSliceRangeZ(&extent[4]);
  inData->GetExtent(inExt);
  int outCheck = 0;
  for (int ii = 0; ii < 3; ii++)
  {
    if (extent[2 * ii] > inExt[2 * ii + 1] || extent[2 * ii + 1] < inExt[2 * ii])
    { // extents don't intersect, we're done
      return;
    }
    if (extent[2 * ii] < inExt[2 * ii])
    {
      extent[2 * ii] = inExt[2 * ii];
    }
    if (extent[2 * ii + 1] > inExt[2 * ii + 1])
    {
      extent[2 * ii + 1] = inExt[2 * ii + 1];
    }
    // check against output extent
    if (extent[2 * ii] < outExt[2 * ii] || extent[2 * ii + 1] > outExt[2 * ii + 1])
    {
      outCheck = 1;
    }
  }

  // Indexing goes from 0 to maxIdX
  int maxIdX = extent[1] - extent[0];
  int maxIdY = extent[3] - extent[2];
  int maxIdZ = extent[5] - extent[4];

  // Convert output limits
  int minOutIdX = outExt[0] - extent[0];
  int maxOutIdX = outExt[1] - extent[0];
  int minOutIdY = outExt[2] - extent[2];
  int maxOutIdY = outExt[3] - extent[2];
  int minOutIdZ = outExt[4] - extent[4];
  int maxOutIdZ = outExt[5] - extent[4];

  // Total number of voxels
  vtkIdType fullsize = (static_cast<vtkIdType>(extent[1] - extent[0] + 1) *
    static_cast<vtkIdType>(extent[3] - extent[2] + 1) *
    static_cast<vtkIdType>(extent[5] - extent[4] + 1));

  // for the progress meter
  double progress = 0.0;
  vtkIdType target = static_cast<vtkIdType>(fullsize / 50.0);
  target++;

  // Setup the mask
  maskData->SetOrigin(inData->GetOrigin());
  maskData->SetSpacing(inData->GetSpacing());
  maskData->SetExtent(extent);
  maskData->AllocateScalars(VTK_UNSIGNED_CHAR, 1);

  unsigned char* maskPtr = static_cast<unsigned char*>(maskData->GetScalarPointerForExtent(extent));
  vtkIdType maskInc[3];
  maskInc[0] = 1;
  maskInc[1] = (extent[1] - extent[0] + 1);
  maskInc[2] = maskInc[1] * (extent[3] - extent[2] + 1);

  // Get input pointer for the extent used by the maskData
  inPtr = static_cast<IT*>(inData->GetScalarPointerForExtent(extent));
  vtkIdType inInc[3];
  inData->GetIncrements(inInc);

  // Get output pointer for the whole output extent
  outPtr = static_cast<OT*>(outData->GetScalarPointerForExtent(outExt));
  vtkIdType outInc[3];
  outData->GetIncrements(outInc);
  // Adjust it so that it corresponds to the maskData extent
  outPtr -= minOutIdX * outInc[0] + minOutIdY * outInc[1] + minOutIdZ * outInc[2];

  // Adjust pointers to active component
  inPtr += activeComponent;
  outPtr += activeComponent;

  if (stencil == nullptr)
  {
    memset(maskPtr, 0, fullsize);
  }
  else
  {
    // pre-set all mask voxels that are outside the stencil
    vtkImageThresholdConnectivityApplyStencil(maskData, stencil, extent);
  }

  // Check whether neighborhood will be used
  double f = self->GetNeighborhoodFraction();
  double radius[3];
  self->GetNeighborhoodRadius(radius);
  int xradius = static_cast<int>(radius[0] + 0.5);
  int yradius = static_cast<int>(radius[1] + 0.5);
  int zradius = static_cast<int>(radius[2] + 0.5);
  double fx = 0.0, fy = 0.0, fz = 0.0;
  bool useNeighborhood = ((xradius > 0) & (yradius > 0) & (zradius > 0));
  if (useNeighborhood)
  {
    fx = 1.0 / radius[0];
    fy = 1.0 / radius[1];
    fz = 1.0 / radius[2];
  }

  // Perform the flood fill within the extent
  double spacing[3];
  double origin[3];
  outData->GetSpacing(spacing);
  outData->GetOrigin(origin);

  // create the seed stack:
  // stack has methods empty(), top(), push(), and pop()
  std::stack<vtkFloodFillSeed> seedStack;

  // initialize with the seeds provided by the user
  vtkPoints* points = self->GetSeedPoints();
  if (points == nullptr)
  { // no seeds!
    return;
  }

  double point[3];
  vtkIdType nPoints = points->GetNumberOfPoints();
  for (vtkIdType p = 0; p < nPoints; p++)
  {
    points->GetPoint(p, point);
    vtkFloodFillSeed seed =
      vtkFloodFillSeed(vtkMath::Floor((point[0] - origin[0]) / spacing[0] + 0.5) - extent[0],
        vtkMath::Floor((point[1] - origin[1]) / spacing[1] + 0.5) - extent[2],
        vtkMath::Floor((point[2] - origin[2]) / spacing[2] + 0.5) - extent[4]);

    if (seed[0] >= 0 && seed[0] <= maxIdX && seed[1] >= 0 && seed[1] <= maxIdY && seed[2] >= 0 &&
      seed[2] <= maxIdZ)
    {
      seedStack.push(seed);
    }
  }

  vtkIdType counter = 0;
  vtkIdType fullcount = 0;

  while (!seedStack.empty())
  {
    vtkFloodFillSeed seed = seedStack.top();
    seedStack.pop();

    unsigned char* maskPtr1 =
      maskPtr + (seed[0] * maskInc[0] + seed[1] * maskInc[1] + seed[2] * maskInc[2]);

    if (*maskPtr1)
    {
      continue;
    }
    *maskPtr1 = 255;

    fullcount++;
    if (id == 0 && (fullcount % target) == 0)
    {
      double v = counter / (10.0 * fullcount);
      double p = fullcount / (v * fullsize + (1.0 - v) * fullcount);
      if (p > progress)
      {
        progress = p;
        self->UpdateProgress(progress);
      }
    }

    IT* inPtr1 = inPtr + (seed[0] * inInc[0] + seed[1] * inInc[1] + seed[2] * inInc[2]);
    IT temp = *inPtr1;

    bool inside = ((lowerThreshold <= temp) & (temp <= upperThreshold));

    // use a spherical neighborhood
    if (useNeighborhood)
    {
      int xmin = seed[0] - xradius;
      xmin = (xmin >= 0 ? xmin : 0);
      int xmax = seed[0] + xradius;
      xmax = (xmax <= maxIdX ? xmax : maxIdX);

      int ymin = seed[1] - yradius;
      ymin = (ymin >= 0 ? ymin : 0);
      int ymax = seed[1] + yradius;
      ymax = (ymax <= maxIdY ? ymax : maxIdY);

      int zmin = seed[2] - zradius;
      zmin = (zmin >= 0 ? zmin : 0);
      int zmax = seed[2] + zradius;
      zmax = (zmax <= maxIdZ ? zmax : maxIdZ);

      inPtr1 = inPtr + (xmin * inInc[0] + ymin * inInc[1] + zmin * inInc[2]);

      int totalcount = 0;
      int threshcount = 0;
      int iz = zmin;
      do
      {
        IT* inPtr2 = inPtr1;
        double rz = (iz - seed[2]) * fz;
        rz *= rz;
        int iy = ymin;
        do
        {
          IT* inPtr3 = inPtr2;
          double ry = (iy - seed[1]) * fy;
          ry *= ry;
          double rzy = rz + ry;
          int ix = xmin;
          do
          {
            double rx = (ix - seed[0]) * fx;
            rx *= rx;
            double rzyx = rzy + rx;
            // include a tolerance in radius check
            bool isgood = (rzyx < (1.0 + 7.62939453125e-06));
            totalcount += isgood;
            isgood &= ((lowerThreshold <= *inPtr3) & (*inPtr3 <= upperThreshold));
            threshcount += isgood;
            inPtr3 += inInc[0];
          } while (++ix <= xmax);
          inPtr2 += inInc[1];
        } while (++iy <= ymax);
        inPtr1 += inInc[2];
      } while (++iz <= zmax);

      // what fraction of the sphere is within threshold?
      inside &= !(static_cast<double>(threshcount) < totalcount * f);
    }

    if (inside)
    {
      // match
      OT* outPtr1 = outPtr + (seed[0] * outInc[0] + seed[1] * outInc[1] + seed[2] * outInc[2]);

      if (outCheck == 0 ||
        (seed[0] >= minOutIdX && seed[0] <= maxOutIdX && seed[1] >= minOutIdY &&
          seed[1] <= maxOutIdY && seed[2] >= minOutIdZ && seed[2] <= maxOutIdZ))
      {
        *outPtr1 = (replaceIn ? inValue : static_cast<OT>(temp));
      }

      // count the seed
      counter += 1;

      // push the new seeds
      if (seed[2] > 0 && *(maskPtr1 - maskInc[2]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0], seed[1], seed[2] - 1));
      }
      if (seed[2] < maxIdZ && *(maskPtr1 + maskInc[2]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0], seed[1], seed[2] + 1));
      }
      if (seed[1] > 0 && *(maskPtr1 - maskInc[1]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0], seed[1] - 1, seed[2]));
      }
      if (seed[1] < maxIdY && *(maskPtr1 + maskInc[1]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0], seed[1] + 1, seed[2]));
      }
      if (seed[0] > 0 && *(maskPtr1 - maskInc[0]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0] - 1, seed[1], seed[2]));
      }
      if (seed[0] < maxIdX && *(maskPtr1 + maskInc[0]) == 0)
      {
        seedStack.push(vtkFloodFillSeed(seed[0] + 1, seed[1], seed[2]));
      }
    }
  }

  if (id == 0)
  {
    self->UpdateProgress(1.0);
  }

  voxelCount = counter;
}

//------------------------------------------------------------------------------
int vtkImageThresholdConnectivity::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* vtkNotUsed(outputVector))
{
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* stencilInfo = inputVector[1]->GetInformationObject(0);

  int inExt[6], extent[6];
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inExt);
  extent[0] = this->SliceRangeX[0];
  extent[1] = this->SliceRangeX[1];
  extent[2] = this->SliceRangeY[0];
  extent[3] = this->SliceRangeY[1];
  extent[4] = this->SliceRangeZ[0];
  extent[5] = this->SliceRangeZ[1];

  // Clip the extent to the inExt
  for (int i = 0; i < 3; i++)
  {
    if (extent[2 * i] < inExt[2 * i])
    {
      extent[2 * i] = inExt[2 * i];
    }
    if (extent[2 * i + 1] > inExt[2 * i + 1])
    {
      extent[2 * i + 1] = inExt[2 * i + 1];
    }
  }

  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent, 6);
  if (stencilInfo)
  {
    stencilInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent, 6);
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkImageThresholdConnectivity::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* stencilInfo = inputVector[1]->GetInformationObject(0);

  vtkImageData* outData = static_cast<vtkImageData*>(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkImageData* inData = static_cast<vtkImageData*>(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkImageData* maskData = this->ImageMask;

  vtkImageStencilData* stencil = nullptr;
  if (stencilInfo)
  {
    stencil = static_cast<vtkImageStencilData*>(stencilInfo->Get(vtkDataObject::DATA_OBJECT()));
  }

  int outExt[6];
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), outExt);
  this->AllocateOutputData(outData, outInfo, outExt);

  // get scalar pointers
  void* inPtr = inData->GetScalarPointerForExtent(outExt);
  void* outPtr = outData->GetScalarPointerForExtent(outExt);

  int id = 0; // not multi-threaded

  if (inData->GetScalarType() != outData->GetScalarType())
  {
    vtkErrorMacro("Execute: Output ScalarType "
      << outData->GetScalarType() << ", must Input ScalarType " << inData->GetScalarType());
    return 0;
  }

  switch (inData->GetScalarType())
  {
    vtkTemplateAliasMacro(
      vtkImageThresholdConnectivityExecute(this, inData, outData, stencil, maskData, outExt, id,
        static_cast<VTK_TT*>(inPtr), static_cast<VTK_TT*>(outPtr), this->NumberOfInVoxels));

    default:
      vtkErrorMacro(<< "Execute: Unknown input ScalarType");
      return 0;
  }

  return 1;
}

//------------------------------------------------------------------------------
void vtkImageThresholdConnectivity::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "InValue: " << this->InValue << "\n";
  os << indent << "OutValue: " << this->OutValue << "\n";
  os << indent << "LowerThreshold: " << this->LowerThreshold << "\n";
  os << indent << "UpperThreshold: " << this->UpperThreshold << "\n";
  os << indent << "ReplaceIn: " << this->ReplaceIn << "\n";
  os << indent << "ReplaceOut: " << this->ReplaceOut << "\n";
  os << indent << "NeighborhoodRadius: " << this->NeighborhoodRadius[0] << " "
     << this->NeighborhoodRadius[1] << " " << this->NeighborhoodRadius[2] << "\n";
  os << indent << "NeighborhoodFraction: " << this->NeighborhoodFraction << "\n";
  os << indent << "NumberOfInVoxels: " << this->NumberOfInVoxels << "\n";
  os << indent << "SliceRangeX: " << this->SliceRangeX[0] << " " << this->SliceRangeX[1] << "\n";
  os << indent << "SliceRangeY: " << this->SliceRangeY[0] << " " << this->SliceRangeY[1] << "\n";
  os << indent << "SliceRangeZ: " << this->SliceRangeZ[0] << " " << this->SliceRangeZ[1] << "\n";
  os << indent << "SeedPoints: " << this->SeedPoints << "\n";
  if (this->SeedPoints)
  {
    this->SeedPoints->PrintSelf(os, indent.GetNextIndent());
  }
  os << indent << "Stencil: " << this->GetStencil() << "\n";
  os << indent << "ActiveComponent: " << this->ActiveComponent << "\n";
}