// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkExtractCTHPart.h"

#include "vtkAppendFilter.h"
#include "vtkAppendPolyData.h"
#include "vtkBoundingBox.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCharArray.h"
#include "vtkClipDataSet.h"
#include "vtkClipPolyData.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkCompositeDataSet.h"
#include "vtkContourFilter.h"
#include "vtkCutter.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkDoubleArray.h"
#include "vtkExecutive.h"
#include "vtkGarbageCollector.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkMultiProcessController.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPlaneCollection.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkRectilinearGrid.h"
#include "vtkSmartPointer.h"
#include "vtkTimerLog.h"
#include "vtkUniformGrid.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <string>
#include <vector>

#include <iostream>

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkExtractCTHPart);
vtkCxxSetObjectMacro(vtkExtractCTHPart, ClipPlane, vtkPlane);
vtkCxxSetObjectMacro(vtkExtractCTHPart, Controller, vtkMultiProcessController);

constexpr double CTH_AMR_SURFACE_VALUE = 0.499;
constexpr double CTH_AMR_SURFACE_VALUE_FLOAT = 1;
constexpr double CTH_AMR_SURFACE_VALUE_UNSIGNED_CHAR = 255;

//------------------------------------------------------------------------------
//=============================================================================
class vtkExtractCTHPartInternal
{
public:
  typedef std::vector<std::string> VolumeArrayNamesType;
  VolumeArrayNamesType VolumeArrayNames;
  vtkBoundingBox GlobalInputBounds;

  // Counter used to scale progress events.
  int TotalNumberOfDatasets;
};

class vtkExtractCTHPart::VectorOfFragments : public std::vector<vtkSmartPointer<vtkPolyData>>
{
};

class vtkExtractCTHPart::VectorOfSolids : public std::vector<vtkSmartPointer<vtkUnstructuredGrid>>
{
};

class vtkExtractCTHPart::ScaledProgress
{
  vtkExtractCTHPart* Self;
  double Shift;
  double Scale;

public:
  ScaledProgress(double shift, double scale, vtkExtractCTHPart* self)
  {
    assert(
      (self != nullptr) && (shift >= 0.0) && (shift <= 1.0) && (scale >= 0.0) && (scale <= 1.0));

    this->Self = self;
    this->Shift = self->ProgressShift;
    this->Scale = self->ProgressScale;

    self->ProgressShift += shift * self->ProgressScale;
    self->ProgressScale *= scale;
    // std::cout << "Shift-Scale Push: " << self->ProgressShift << ", " <<
    //  self->ProgressScale << endl;
  }

  ~ScaledProgress() { this->WorkDone(); }

  void WorkDone()
  {
    if (this->Self)
    {
      this->Self->ProgressScale = this->Scale;
      this->Self->ProgressShift = this->Shift;
      // std::cout << "Shift-Scale Pop: " << this->Self->ProgressShift << ", " <<
      //  this->Self->ProgressScale << endl;
      this->Self = nullptr;
    }
  }
};

//=============================================================================
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
vtkExtractCTHPart::vtkExtractCTHPart()
{
  this->Internals = new vtkExtractCTHPartInternal();
  this->ClipPlane = nullptr;
  this->GenerateTriangles = true;
  this->GenerateSolidGeometry = false;
  this->Capping = true;
  this->RemoveGhostCells = true;
  this->VolumeFractionSurfaceValueInternal = CTH_AMR_SURFACE_VALUE;
  this->VolumeFractionSurfaceValue = CTH_AMR_SURFACE_VALUE;
  this->ProgressScale = 1.0;
  this->ProgressShift = 0.0;

  this->Controller = nullptr;
  this->SetController(vtkMultiProcessController::GetGlobalController());
}

//------------------------------------------------------------------------------
vtkExtractCTHPart::~vtkExtractCTHPart()
{
  this->SetController(nullptr);
  this->SetClipPlane(nullptr);

  delete this->Internals;
  this->Internals = nullptr;
}

//------------------------------------------------------------------------------
// Overload standard modified time function. If clip plane is modified,
// then this object is modified as well.
vtkMTimeType vtkExtractCTHPart::GetMTime()
{
  vtkMTimeType mTime = this->Superclass::GetMTime();
  if (this->ClipPlane)
  {
    vtkMTimeType time = this->ClipPlane->GetMTime();
    return time > mTime ? time : mTime;
  }

  return mTime;
}

//------------------------------------------------------------------------------
void vtkExtractCTHPart::RemoveVolumeArrayNames()
{
  this->Internals->VolumeArrayNames.clear();
  this->Modified();
}

//------------------------------------------------------------------------------
void vtkExtractCTHPart::AddVolumeArrayName(const char* arrayName)
{
  if (arrayName != nullptr && arrayName[0] != 0 &&
    std::find(this->Internals->VolumeArrayNames.begin(), this->Internals->VolumeArrayNames.end(),
      std::string(arrayName)) == this->Internals->VolumeArrayNames.end())
  {
    this->Internals->VolumeArrayNames.emplace_back(arrayName);

    // ensure that the volume arrays are in determinate order. I should just
    // change the code to use a std::set.
    std::sort(this->Internals->VolumeArrayNames.begin(), this->Internals->VolumeArrayNames.end());
    this->Modified();
  }
}

//------------------------------------------------------------------------------
int vtkExtractCTHPart::GetNumberOfVolumeArrayNames()
{
  return static_cast<int>(this->Internals->VolumeArrayNames.size());
}

//------------------------------------------------------------------------------
const char* vtkExtractCTHPart::GetVolumeArrayName(int idx)
{
  if (idx < 0 || idx > static_cast<int>(this->Internals->VolumeArrayNames.size()))
  {
    return nullptr;
  }

  return this->Internals->VolumeArrayNames[idx].c_str();
}

//------------------------------------------------------------------------------
int vtkExtractCTHPart::FillInputPortInformation(int port, vtkInformation* info)
{
  if (!this->Superclass::FillInputPortInformation(port, info))
  {
    return 0;
  }

  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkNonOverlappingAMR");
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkMultiBlockDataSet");
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid");
  return 1;
}

//------------------------------------------------------------------------------
int vtkExtractCTHPart::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  const int number_of_volume_arrays = static_cast<int>(this->Internals->VolumeArrayNames.size());
  if (number_of_volume_arrays == 0)
  {
    // nothing to do.
    return 1;
  }

  vtkDataObject* inputDO = vtkDataObject::GetData(inputVector[0], 0);
  vtkSmartPointer<vtkCompositeDataSet> inputCD = vtkCompositeDataSet::SafeDownCast(inputDO);
  vtkRectilinearGrid* inputRG = vtkRectilinearGrid::SafeDownCast(inputDO);
  assert(inputCD != nullptr || inputRG != nullptr);

  if (inputRG)
  {
    vtkNew<vtkMultiBlockDataSet> mb;
    mb->SetBlock(0, inputRG);
    inputCD = mb;
  }

  vtkMultiBlockDataSet* output = vtkMultiBlockDataSet::GetData(outputVector, 0);

  // initialize output multiblock-dataset. It will always have as many blocks as
  // the number-of-volume arrays requested.
  output->SetNumberOfBlocks(number_of_volume_arrays);

  // Compute global bounds for the input dataset. This is used to generate
  // external surface for the dataset.
  if (!this->ComputeGlobalBounds(inputCD))
  {
    vtkErrorMacro("Failed to compute global bounds.");
    return 0;
  }

  if (!this->Internals->GlobalInputBounds.IsValid())
  {
    // empty input, do nothing.
    return 1;
  }

  unsigned int array_index = 0;
  for (vtkExtractCTHPartInternal::VolumeArrayNamesType::iterator iter =
         this->Internals->VolumeArrayNames.begin();
       iter != this->Internals->VolumeArrayNames.end(); ++iter, ++array_index)
  {
    // this loop is doing the 1/(num-arrays)'th work for the entire task.
    ScaledProgress sp(array_index * 1.0 / this->Internals->VolumeArrayNames.size(),
      1.0 / this->Internals->VolumeArrayNames.size(), this);

    output->GetMetaData(array_index)->Set(vtkCompositeDataSet::NAME(), iter->c_str());

    vtkSmartPointer<vtkDataSet> part = nullptr;
    vtkGarbageCollector::DeferredCollectionPush();

    if (this->GenerateSolidGeometry)
    {
      part = this->ExtractSolid(inputCD, iter->c_str());
    }
    else
    {
      part = this->ExtractContour(inputCD, iter->c_str());
    }

    if (part != nullptr && part->GetNumberOfPoints() > 0)
    {
      // Add extra arrays.
      vtkNew<vtkIntArray> partArray;
      partArray->SetName("Part Index");
      partArray->SetNumberOfComponents(1);
      partArray->SetNumberOfTuples(part->GetNumberOfPoints());
      partArray->FillComponent(0, static_cast<double>(array_index));
      part->GetPointData()->AddArray(partArray);

      // I'm not adding the "Name" array that was added in previous
      // implementation. Don't think that's much of use.

      output->SetBlock(array_index, part);
    }
    vtkGarbageCollector::DeferredCollectionPop();
  }
  return 1;
}

//------------------------------------------------------------------------------
bool vtkExtractCTHPart::ComputeGlobalBounds(vtkCompositeDataSet* input)
{
  assert("pre: input_exists" && input != nullptr);
  this->Internals->GlobalInputBounds.Reset();

  this->Internals->TotalNumberOfDatasets = 0;

  vtkCompositeDataIterator* iter = input->NewIterator();
  for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
  {
    vtkDataSet* ds = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
    if (!ds) // can be null if on another processor
    {
      continue;
    }
    double realBounds[6];
    ds->GetBounds(realBounds);
    this->Internals->GlobalInputBounds.AddBounds(realBounds);

    ++this->Internals->TotalNumberOfDatasets;
  }
  iter->Delete();

  // Here we have the bounds according to our local datasets.
  // If we are not running in parallel then the local
  // bounds are the global bounds
  if (!this->Controller || this->Controller->GetNumberOfProcesses() <= 1)
  {
    return true;
  }

  const double* min_point = this->Internals->GlobalInputBounds.GetMinPoint();
  const double* max_point = this->Internals->GlobalInputBounds.GetMaxPoint();
  double min_result[3], max_result[3];

  if (!this->Controller->AllReduce(min_point, min_result, 3, vtkCommunicator::MIN_OP))
  {
    return false;
  }
  if (!this->Controller->AllReduce(max_point, max_result, 3, vtkCommunicator::MAX_OP))
  {
    return false;
  }

  this->Internals->GlobalInputBounds.SetBounds(
    min_result[0], max_result[0], min_result[1], max_result[1], min_result[2], max_result[2]);
  // At this point, the global bounds is set in each processor.
  return true;
}

//------------------------------------------------------------------------------
// return nullptr on error.
vtkSmartPointer<vtkDataSet> vtkExtractCTHPart::ExtractContour(
  vtkCompositeDataSet* input, const char* arrayName)
{
  assert(input != nullptr && arrayName != nullptr && arrayName[0] != 0);

  bool warn_once = true;
  vtkSmartPointer<vtkCompositeDataIterator> iter;
  iter.TakeReference(input->NewIterator());

  // this loop is first 95% of the work.
  ScaledProgress sp1(0.0, 0.95, this);

  int counter = 0;
  vtkExtractCTHPart::VectorOfFragments fragments;
  for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem(), ++counter)
  {
    // each iteration is 1/(total num of datasets)'th for the work.
    ScaledProgress sp(counter * 1.0 / this->Internals->TotalNumberOfDatasets,
      1.0 / this->Internals->TotalNumberOfDatasets, this);

    if (counter % 1000 == 0)
    {
      this->TriggerProgressEvent(0.0);
    }

    vtkDataObject* dataObj = iter->GetCurrentDataObject();
    vtkRectilinearGrid* rg = vtkRectilinearGrid::SafeDownCast(dataObj);
    vtkUniformGrid* ug = vtkUniformGrid::SafeDownCast(dataObj);

    if (ug)
    {
      if (!this->ExtractClippedContourOnBlock<vtkUniformGrid>(fragments, ug, arrayName))
      {
        return nullptr;
      }
    }
    else if (rg)
    {
      if (!this->ExtractClippedContourOnBlock<vtkRectilinearGrid>(fragments, rg, arrayName))
      {
        return nullptr;
      }
    }
    else if (warn_once && dataObj)
    {
      warn_once = false;
      vtkWarningMacro(<< dataObj->GetClassName() << " will be ignored.");
    }
    if ((counter % 1000) == 0)
    {
      this->TriggerProgressEvent(1.0);
    }
  }

  vtkSmartPointer<vtkDataSet> output = vtkSmartPointer<vtkPolyData>::New();
  if (fragments.empty())
  {
    // empty contour. Not an error though, hence we don't return nullptr.
    return output;
  }
  sp1.WorkDone();

  // Now, the last .05 % of the work.
  ScaledProgress sp2(0.95, 0.05, this);
  this->TriggerProgressEvent(0.0);
  vtkNew<vtkAppendPolyData> appender;
  for (size_t cc = 0; cc < fragments.size(); cc++)
  {
    appender->AddInputData(fragments[cc]);
  }
  appender->Update();

  output->ShallowCopy(appender->GetOutputDataObject(0));
  this->TriggerProgressEvent(1.0);
  return output;
}

//------------------------------------------------------------------------------
// return nullptr on error.
vtkSmartPointer<vtkDataSet> vtkExtractCTHPart::ExtractSolid(
  vtkCompositeDataSet* input, const char* arrayName)
{
  assert(input != nullptr && arrayName != nullptr && arrayName[0] != 0);

  bool warn_once = true;
  vtkSmartPointer<vtkCompositeDataIterator> iter;
  iter.TakeReference(input->NewIterator());

  // this loop is first 95% of the work.
  ScaledProgress sp1(0.0, 0.95, this);

  int counter = 0;
  vtkExtractCTHPart::VectorOfSolids solids;
  for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem(), ++counter)
  {
    // each iteration is 1/(total num of datasets)'th for the work.
    ScaledProgress sp(counter * 1.0 / this->Internals->TotalNumberOfDatasets,
      1.0 / this->Internals->TotalNumberOfDatasets, this);

    if (counter % 1000 == 0)
    {
      this->TriggerProgressEvent(0.0);
    }

    vtkDataObject* dataObj = iter->GetCurrentDataObject();
    vtkRectilinearGrid* rg = vtkRectilinearGrid::SafeDownCast(dataObj);
    vtkUniformGrid* ug = vtkUniformGrid::SafeDownCast(dataObj);

    if (ug)
    {
      if (!this->ExtractClippedVolumeOnBlock<vtkUniformGrid>(solids, ug, arrayName))
      {
        return nullptr;
      }
    }
    else if (rg)
    {
      if (!this->ExtractClippedVolumeOnBlock<vtkRectilinearGrid>(solids, rg, arrayName))
      {
        return nullptr;
      }
    }
    else if (warn_once && dataObj)
    {
      warn_once = false;
      vtkWarningMacro(<< dataObj->GetClassName() << " will be ignored.");
    }
    if ((counter % 1000) == 0)
    {
      this->TriggerProgressEvent(1.0);
    }
  }

  vtkSmartPointer<vtkDataSet> output = vtkSmartPointer<vtkUnstructuredGrid>::New();
  if (solids.empty())
  {
    // empty fragments. Not an error though, hence we don't return nullptr.
    return output;
  }
  sp1.WorkDone();

  // Now, the last .05 % of the work.
  ScaledProgress sp2(0.95, 0.05, this);
  this->TriggerProgressEvent(0.0);
  vtkNew<vtkAppendFilter> appender;
  for (size_t cc = 0; cc < solids.size(); cc++)
  {
    appender->AddInputData(solids[cc]);
  }
  appender->Update();

  output->ShallowCopy(appender->GetOutputDataObject(0));
  this->TriggerProgressEvent(1.0);
  return output;
}

void vtkExtractCTHPart::DetermineSurfaceValue(int dataType)
{
  // determine the true value to use for the contour based on the data-type.
  switch (dataType)
  {
    case VTK_UNSIGNED_CHAR:
      this->VolumeFractionSurfaceValueInternal =
        CTH_AMR_SURFACE_VALUE_UNSIGNED_CHAR * this->VolumeFractionSurfaceValue;
      break;

    default:
      this->VolumeFractionSurfaceValueInternal =
        CTH_AMR_SURFACE_VALUE_FLOAT * this->VolumeFractionSurfaceValue;
  }
}

//------------------------------------------------------------------------------
template <class T>
bool vtkExtractCTHPart::ExtractClippedContourOnBlock(
  vtkExtractCTHPart::VectorOfFragments& fragments, T* dataset, const char* arrayName)
{
  assert(arrayName != nullptr && arrayName[0] != 0 && dataset != nullptr);

  vtkDataArray* volumeFractionArray = dataset->GetCellData()->GetArray(arrayName);
  if (!volumeFractionArray)
  {
    // skip this block.
    return true;
  }

  // determine the true value to use for the contour based on the data-type.
  this->DetermineSurfaceValue(volumeFractionArray->GetDataType());

  // We create a clone so we can modify the dataset (i.e. add new arrays to it).
  vtkNew<T> inputClone;
  inputClone->ShallowCopy(dataset);

  // Convert cell-data-2-point-data so we can contour.
  vtkNew<vtkDoubleArray> pointVolumeFractionArray;
  this->ExecuteCellDataToPointData(
    volumeFractionArray, pointVolumeFractionArray, inputClone->GetDimensions());
  inputClone->GetPointData()->SetScalars(pointVolumeFractionArray);

  VectorOfFragments blockFragments;
  if (!this->ExtractContourOnBlock<T>(blockFragments, inputClone, arrayName))
  {
    return false;
  }

  if (!this->ClipPlane)
  {
    fragments.insert(fragments.end(), blockFragments.begin(), blockFragments.end());
    return true;
  }

  // Clip-n-cap the fragments using the clip plane.

  // for the clip.
  for (size_t cc = 0; cc < blockFragments.size(); cc++)
  {
    vtkNew<vtkClipPolyData> clipper;
    clipper->SetClipFunction(this->ClipPlane);
    clipper->SetInputDataObject(blockFragments[cc]);
    clipper->Update();
    fragments.push_back(clipper->GetOutput());
  }

  //// for the cap.
  if (this->Capping)
  {
    vtkNew<vtkCutter> cutter;
    cutter->SetCutFunction(this->ClipPlane);
    cutter->SetGenerateTriangles(this->GenerateTriangles ? 1 : 0);
    cutter->SetInputDataObject(inputClone);

    vtkNew<vtkClipPolyData> scalarClipper;
    scalarClipper->SetInputConnection(cutter->GetOutputPort());
    scalarClipper->SetValue(this->VolumeFractionSurfaceValueInternal);
    scalarClipper->Update();
    fragments.push_back(scalarClipper->GetOutput());
  }
  return true;
}

//------------------------------------------------------------------------------
template <class T>
bool vtkExtractCTHPart::ExtractContourOnBlock(
  vtkExtractCTHPart::VectorOfFragments& fragments, T* dataset, const char* arrayName)
{
  assert(arrayName != nullptr && arrayName[0] != 0 && dataset != nullptr);

  vtkPointData* pd = dataset->GetPointData();

  // Contour only if necessary.
  double range[2];
  if (!pd->GetRange(arrayName, range))
  {
    return false;
  }
  if (range[1] < this->VolumeFractionSurfaceValueInternal)
  {
    // this block doesn't have the material of interest.
    return true;
  }

  // Extract exterior surface. Adds the surface polydata to fragments, if any.
  if (this->Capping)
  {
    this->ExtractExteriorSurface(fragments, dataset);
  }

  if (this->ClipPlane == nullptr && range[0] > this->VolumeFractionSurfaceValueInternal)
  {
    // no need to extract contour.
    return true;
  }

  // Extract contour.
  vtkNew<vtkContourFilter> contourer;
  contourer->SetInputData(dataset);
  contourer->SetValue(0, this->VolumeFractionSurfaceValueInternal);
  contourer->SetComputeScalars(0);
  contourer->SetGenerateTriangles(this->GenerateTriangles ? 1 : 0);
  contourer->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, arrayName);
  contourer->Update();

  vtkPolyData* output = vtkPolyData::SafeDownCast(contourer->GetOutputDataObject(0));
  if (!output || output->GetNumberOfPoints() == 0)
  {
    return true;
  }
  if (!this->RemoveGhostCells)
  {
    // BUG #14291. Rather than renaming the array, we remove the GhostArray
    // from the output since it may not be present on all ranks and will cause array count mismatch
    output->GetCellData()->RemoveArray(vtkDataSetAttributes::GhostArrayName());
  }

  fragments.push_back(output);
  return true;
}

//------------------------------------------------------------------------------
// Description:
// Append quads for faces of the block that actually on the bounds
// of the hierarchical dataset. Deals with ghost cells.
template <class T>
void vtkExtractCTHPart::ExtractExteriorSurface(
  vtkExtractCTHPart::VectorOfFragments& fragments, T* input)
{
  assert("pre: valid_input" && input != nullptr && input->CheckAttributes() == 0);

  int result = 0;
#if 1
  int dims[3];
  input->GetDimensions(dims);
  int ext[6];
  int originalExtents[6];
  input->GetExtent(ext);
  input->GetExtent(originalExtents);

  //  vtkUnsignedCharArray *ghostArray=static_cast<vtkUnsignedCharArray
  //  *>(input->GetCellData()->GetArray(vtkDataSetAttributes::GhostArrayName()));

  // bounds without taking ghost cells into account
  double bounds[6];

  input->GetBounds(bounds);

#if 0
  // block face min x
  if(this->IsGhostFace(0,0,dims,ghostArray))
  {
    // downsize this!
    bounds[0]+=spacing[0];
    ++ext[0];
  }
  if(this->IsGhostFace(0,1,dims,ghostArray))
  {
    // downsize this!
    bounds[1]-=spacing[0];
    --ext[1];
  }
  if(this->IsGhostFace(1,0,dims,ghostArray))
  {
    // downsize this!
    bounds[2]+=spacing[1];
    ++ext[2];
  }
  if(this->IsGhostFace(1,1,dims,ghostArray))
  {
    // downsize this!
    bounds[3]-=spacing[1];
    --ext[3];
  }
  if(this->IsGhostFace(2,0,dims,ghostArray))
  {
    // downsize this!
    bounds[4]+=spacing[2];
    ++ext[4];
  }
  if(this->IsGhostFace(2,1,dims,ghostArray))
  {
    // downsize this!
    bounds[5]-=spacing[2];
    --ext[5];
  }
#endif
  // here, bounds are real block bounds without ghostcells.

  const double* minP = this->Internals->GlobalInputBounds.GetMinPoint();
  const double* maxP = this->Internals->GlobalInputBounds.GetMaxPoint();
#if 0
  const double epsilon=0.001;
  int doFaceMinX=fabs(bounds[0]- minP[0])<epsilon;
  int doFaceMaxX=fabs(bounds[1]- maxP[0])<epsilon;
  int doFaceMinY=fabs(bounds[2]- minP[1])<epsilon;
  int doFaceMaxY=fabs(bounds[3]- maxP[1])<epsilon;
  int doFaceMinZ=fabs(bounds[4]- minP[2])<epsilon;
  int doFaceMaxZ=fabs(bounds[5]- maxP[2])<epsilon;
#endif

#if 1
  int doFaceMinX = bounds[0] <= minP[0];
  int doFaceMaxX = bounds[1] >= maxP[0];
  int doFaceMinY = bounds[2] <= minP[1];
  int doFaceMaxY = bounds[3] >= maxP[1];
  int doFaceMinZ = bounds[4] <= minP[2];
  int doFaceMaxZ = bounds[5] >= maxP[2];
#endif
#if 0
  int doFaceMinX=1;
  int doFaceMaxX=1;
  int doFaceMinY=1;
  int doFaceMaxY=1;
  int doFaceMinZ=1;
  int doFaceMaxZ=1;
#endif
#if 0
  int doFaceMinX=0;
  int doFaceMaxX=0;
  int doFaceMinY=0;
  int doFaceMaxY=0;
  int doFaceMinZ=0;
  int doFaceMaxZ=0;
#endif
#if 0
  doFaceMaxX=0;
  doFaceMaxY=0;
  doFaceMaxZ=0;
#endif

  result = doFaceMinX || doFaceMaxX || doFaceMinY || doFaceMaxY || doFaceMinZ || doFaceMaxZ;

  if (result)
  {
    vtkSmartPointer<vtkPolyData> output = vtkSmartPointer<vtkPolyData>::New();

    vtkIdType numPoints = 0;
    vtkIdType cellArraySize = 0;

    //  input->GetExtent(ext);

    // Compute an upper bound for the number of points and cells.
    // xMin face
    if (doFaceMinX && ext[2] != ext[3] && ext[4] != ext[5] && ext[0] != ext[1])
    {
      cellArraySize += 2 * (ext[3] - ext[2] + 1) * (ext[5] - ext[4] + 1);
      numPoints += (ext[3] - ext[2] + 1) * (ext[5] - ext[4] + 1);
    }
    // xMax face
    if (doFaceMaxX && ext[2] != ext[3] && ext[4] != ext[5])
    {
      cellArraySize += 2 * (ext[3] - ext[2] + 1) * (ext[5] - ext[4] + 1);
      numPoints += (ext[3] - ext[2] + 1) * (ext[5] - ext[4] + 1);
    }
    // yMin face
    if (doFaceMinY && ext[0] != ext[1] && ext[4] != ext[5] && ext[2] != ext[3])
    {
      cellArraySize += 2 * (ext[1] - ext[0] + 1) * (ext[5] - ext[4] + 1);
      numPoints += (ext[1] - ext[0] + 1) * (ext[5] - ext[4] + 1);
    }
    // yMax face
    if (doFaceMaxY && ext[0] != ext[1] && ext[4] != ext[5])
    {
      cellArraySize += 2 * (ext[1] - ext[0] + 1) * (ext[5] - ext[4] + 1);
      numPoints += (ext[1] - ext[0] + 1) * (ext[5] - ext[4] + 1);
    }
    // zMin face
    if (doFaceMinZ && ext[0] != ext[1] && ext[2] != ext[3] && ext[4] != ext[5])
    {
      cellArraySize += 2 * (ext[1] - ext[0] + 1) * (ext[3] - ext[2] + 1);
      numPoints += (ext[1] - ext[0] + 1) * (ext[3] - ext[2] + 1);
    }
    // zMax face
    if (doFaceMaxZ && ext[0] != ext[1] && ext[2] != ext[3])
    {
      cellArraySize += 2 * (ext[1] - ext[0] + 1) * (ext[3] - ext[2] + 1);
      numPoints += (ext[1] - ext[0] + 1) * (ext[3] - ext[2] + 1);
    }

    vtkCellArray* outPolys = vtkCellArray::New();
    outPolys->AllocateEstimate(cellArraySize, 1);
    output->SetPolys(outPolys);
    outPolys->Delete();

    vtkPoints* outPoints = vtkPoints::New();
    outPoints->Allocate(numPoints);
    output->SetPoints(outPoints);
    outPoints->Delete();

    // Allocate attributes for copying.
    output->GetPointData()->CopyAllocate(input->GetPointData());
    output->GetCellData()->CopyAllocate(input->GetCellData());

    // Extents are already corrected for ghostcells.

    // make each face that is actually on the ds boundary
    if (doFaceMinX)
    {
      this->ExecuteFaceQuads(input, output, 0, originalExtents, ext, 0, 1, 2);
    }
    if (doFaceMaxX)
    {
      this->ExecuteFaceQuads(input, output, 1, originalExtents, ext, 0, 2, 1);
    }
    if (doFaceMinY)
    {
      this->ExecuteFaceQuads(input, output, 0, originalExtents, ext, 1, 2, 0);
    }
    if (doFaceMaxY)
    {
      this->ExecuteFaceQuads(input, output, 1, originalExtents, ext, 1, 0, 2);
    }
    if (doFaceMinZ)
    {
      this->ExecuteFaceQuads(input, output, 0, originalExtents, ext, 2, 0, 1);
    }
    if (doFaceMaxZ)
    {
      this->ExecuteFaceQuads(input, output, 1, originalExtents, ext, 2, 1, 0);
    }
    output->Squeeze();
    assert(output->CheckAttributes() == 0);

    vtkNew<vtkClipPolyData> clipper;
    clipper->SetInputData(output);
    clipper->SetValue(this->VolumeFractionSurfaceValueInternal);
    clipper->Update();
    fragments.push_back(clipper->GetOutput());
  }
#endif
  // result=>valid_surface: A=>B !A||B
}

template <class T>
bool vtkExtractCTHPart::ExtractClippedVolumeOnBlock(
  VectorOfSolids& solids, T* dataset, const char* arrayName)
{
  assert(arrayName != nullptr && arrayName[0] != 0 && dataset != nullptr);

  vtkCellData* cd = dataset->GetCellData();
  vtkDataArray* volumeFractionArray = cd->GetArray(arrayName);
  if (!volumeFractionArray)
  {
    // skip this block.
    return true;
  }

  // determine the true value to use for clipping based on the data-type.
  this->DetermineSurfaceValue(volumeFractionArray->GetDataType());

  // We create a clone so we can modify the dataset (i.e. add new arrays to it).
  vtkNew<T> inputClone;
  inputClone->ShallowCopy(dataset);

  // Convert cell-data-2-point-data so we can clip the cells.
  vtkNew<vtkDoubleArray> pointVolumeFractionArray;
  this->ExecuteCellDataToPointData(
    volumeFractionArray, pointVolumeFractionArray, inputClone->GetDimensions());
  inputClone->GetPointData()->SetScalars(pointVolumeFractionArray);

  // clip volume only if necessary.
  double range[2];
  if (!cd->GetRange(arrayName, range))
  {
    return false;
  }
  if (range[0] > this->VolumeFractionSurfaceValueInternal ||
    range[1] < this->VolumeFractionSurfaceValueInternal)
  {
    // this block doesn't have the material of interest.
    return true;
  }

  // Clip the volume.
  vtkNew<vtkClipDataSet> blockClipper;
  blockClipper->SetInputData(inputClone);
  blockClipper->SetValue(this->VolumeFractionSurfaceValueInternal);
  blockClipper->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, arrayName);
  blockClipper->Update();
  vtkSmartPointer<vtkUnstructuredGrid> solidFragment = blockClipper->GetOutput();

  solidFragment->GetPointData()->RemoveArray(arrayName);

  if (!this->ClipPlane)
  {
    solids.push_back(solidFragment);
    return true;
  }

  // Clip the solid fragment using the user-specified clip plane.
  vtkNew<vtkClipDataSet> fragClipper;
  fragClipper->SetClipFunction(this->ClipPlane);
  fragClipper->SetInputData(solidFragment);
  fragClipper->Update();
  solids.push_back(fragClipper->GetOutput());

  return true;
}

//------------------------------------------------------------------------------
// Description:
// Is block face on axis0 (either min or max depending on the maxFlag)
// composed of only ghost cells?
int vtkExtractCTHPart::IsGhostFace(
  int axis0, int maxFlag, int dims[3], vtkUnsignedCharArray* ghostArray)
{
  assert("pre: valid_axis0" && axis0 >= 0 && axis0 <= 2);

  int axis1 = axis0 + 1;
  if (axis1 > 2)
  {
    axis1 = 0;
  }
  int axis2 = axis0 + 2;
  if (axis2 > 2)
  {
    axis2 = 0;
  }

  int ijk[3]; // index of the cell.

  if (maxFlag)
  {
    ijk[axis0] = dims[axis0] - 2;
  }
  else
  {
    ijk[axis0] = 0;
  }

  // We test the center cell of the block face.
  // in the worst case (2x2 cells), we need to know if at least
  // three cells are ghost to say that the face is a ghost face.

  ijk[axis1] = dims[axis1] / 2 - 1; // (dims[axis1]-2)/2
  ijk[axis2] = dims[axis2] / 2 - 1; // (dims[axis2]-2)/2
  int result = ghostArray->GetValue(vtkStructuredData::ComputeCellId(dims, ijk));

  if (dims[axis1] == 3)
  {
    // axis1 requires 2 cells to be tested.
    // if so, axis1index=0 and axis1index+1=1
    ijk[axis1] = 1;
    result = result && ghostArray->GetValue(vtkStructuredData::ComputeCellId(dims, ijk));
  }

  if (dims[axis2] == 3)
  {
    // herem axis1 may have moved from the previous test.
    // axis2 requires 2 cells to be tested.
    // if so, axis2index=0 and axis2index+1=1
    ijk[axis2] = 1;
    result = result && ghostArray->GetValue(vtkStructuredData::ComputeCellId(dims, ijk));
  }
  return result;
}

//------------------------------------------------------------------------------
// Description:
// Merely the same implementation than in vtkDataSetSurfaceFilter, without
// dealing with the whole extents.
void vtkExtractCTHPart::ExecuteFaceQuads(vtkDataSet* input, vtkPolyData* output, int maxFlag,
  int originalExtents[6], int ext[6], int aAxis, int bAxis, int cAxis)
{
  assert("pre: input_exists" && input != nullptr);
  assert("pre: output_exists" && output != nullptr);
  assert("pre: originalExtents_exists" && originalExtents != nullptr);
  assert("pre: ext_exists" && ext != nullptr);
  assert("pre: valid_axes" && aAxis >= 0 && aAxis <= 2 && bAxis >= 0 && bAxis <= 2 && cAxis >= 0 &&
    cAxis <= 2 && aAxis != bAxis && aAxis != cAxis && bAxis != cAxis);

  vtkPoints* outPts;
  vtkCellArray* outPolys;
  vtkPointData *inPD, *outPD;
  vtkCellData *inCD, *outCD;
  int pInc[3];
  int qInc[3];
  int cOutInc;
  double pt[3];
  vtkIdType inStartPtId;
  vtkIdType inStartCellId;
  vtkIdType outStartPtId;
  vtkIdType outPtId;
  vtkIdType inId, outId;
  int ib, ic;
  int aA2, bA2, cA2;

  outPts = output->GetPoints();
  outPD = output->GetPointData();
  inPD = input->GetPointData();
  outCD = output->GetCellData();
  inCD = input->GetCellData();

  pInc[0] = 1;
  pInc[1] = (originalExtents[1] - originalExtents[0] + 1);
  pInc[2] = (originalExtents[3] - originalExtents[2] + 1) * pInc[1];
  // quad increments (cell increments, but cInc could be confused with c axis).
  qInc[0] = 1;
  qInc[1] = originalExtents[1] - originalExtents[0];
  // The conditions are for when we have one or more degenerate axes (2d or 1d cells).
  if (qInc[1] == 0)
  {
    qInc[1] = 1;
  }
  qInc[2] = (originalExtents[3] - originalExtents[2]) * qInc[1];
  if (qInc[2] == 0)
  {
    qInc[2] = qInc[1];
  }

  // Temporary variables to avoid many multiplications.
  aA2 = aAxis << 1; // * 2;
  bA2 = bAxis << 1; // * 2;
  cA2 = cAxis << 1; //  * 2;

  // We might as well put the test for this face here.
  if (ext[bA2] == ext[bA2 + 1] || ext[cA2] == ext[cA2 + 1])
  {
    return;
  }
#if 0
  if (maxFlag)
  {
    if (ext[aA2+1] < wholeExt[aA2+1])
    {
      return;
    }
  }
  else
  { // min faces have a slightly different condition to avoid coincident faces.
    if (ext[aA2] == ext[aA2+1] || ext[aA2] > wholeExt[aA2])
    {
      return;
    }
  }
#endif

  if (!maxFlag)
  {
    if (ext[aA2] == ext[aA2 + 1])
    {
      return;
    }
  }

  // Assuming no ghost cells ...
  //  inStartPtId = inStartCellId = 0;
  inStartPtId = 0;   // ext[aA2];
  inStartCellId = 0; // ext[aA2];

  // I put this confusing conditional to fix a regression test.
  // If we are creating a maximum face, then we indeed have to offset the input cell Ids.
  // However, vtkGeometryFilter created a 2d image as a max face, but the cells are copied
  // as a min face (no offset).  Hence maxFlag = 1 and there should be no offset.
  if (maxFlag && ext[aA2] < ext[1 + aA2])
  {
    inStartPtId = pInc[aAxis] * (ext[aA2 + 1] - originalExtents[aA2]);       // -ext[aA2]
    inStartCellId = qInc[aAxis] * (ext[aA2 + 1] - originalExtents[aA2] - 1); // -ext[aA2]
  }

  outStartPtId = outPts->GetNumberOfPoints();
  // Make the points for this face.
  for (ic = ext[cA2]; ic <= ext[cA2 + 1]; ++ic)
  {
    for (ib = ext[bA2]; ib <= ext[bA2 + 1]; ++ib)
    {
      //      inId = inStartPtId + (ib-ext[bA2]+originExtents[bAxis])*pInc[bAxis]
      //                         + (ic-ext[cA2]+originExtents[cAxis])*pInc[cAxis];

      inId = inStartPtId + (ib - originalExtents[bA2]) * pInc[bAxis] +
        (ic - originalExtents[cA2]) * pInc[cAxis];

      input->GetPoint(inId, pt);
      outId = outPts->InsertNextPoint(pt);
      // Copy point data.
      outPD->CopyData(inPD, inId, outId);
    }
  }

  // Do the cells.
  cOutInc = ext[bA2 + 1] - ext[bA2] + 1;

  outPolys = output->GetPolys();

  // Old method for creating quads (needed for cell data.).
  for (ic = ext[cA2]; ic < ext[cA2 + 1]; ++ic)
  {
    for (ib = ext[bA2]; ib < ext[bA2 + 1]; ++ib)
    {
      outPtId = outStartPtId + (ib - ext[bA2]) + (ic - ext[cA2]) * cOutInc;
      //      inId = inStartCellId + (ib-ext[bA2]+originExtents[bAxis])*qInc[bAxis] +
      //      (ic-ext[cA2]+originExtents[cAxis])*qInc[cAxis];

      inId = inStartCellId + (ib - originalExtents[bA2]) * qInc[bAxis] +
        (ic - originalExtents[cA2]) * qInc[cAxis];

      outId = outPolys->InsertNextCell(4);
      outPolys->InsertCellPoint(outPtId);
      outPolys->InsertCellPoint(outPtId + cOutInc);
      outPolys->InsertCellPoint(outPtId + cOutInc + 1);
      outPolys->InsertCellPoint(outPtId + 1);

      // Copy cell data.
      outCD->CopyData(inCD, inId, outId);
    }
  }
}

//------------------------------------------------------------------------------
void vtkExtractCTHPart::ExecuteCellDataToPointData(
  vtkDataArray* cellVolumeFraction, vtkDoubleArray* pointVolumeFraction, const int* dims)
{
  int count;
  int i, j, k;
  int iEnd, jEnd, kEnd;
  int jInc, kInc;
  double* pPoint;

  pointVolumeFraction->SetName(cellVolumeFraction->GetName());
  pointVolumeFraction->SetNumberOfTuples(dims[0] * dims[1] * dims[2]);

  iEnd = dims[0] - 1;
  jEnd = dims[1] - 1;
  kEnd = dims[2] - 1;

  // Deals with none 3D images, otherwise it will never enter into the loop.
  // And then the data will be not initialized and the output of the contour
  // will be empty.

  int dimensionality = 3;

  if (kEnd == 0)
  {
    --dimensionality;
    kEnd = 1;
  }

  // Increments are for the point array.
  jInc = dims[0];
  kInc = (dims[1]) * jInc;

  pPoint = pointVolumeFraction->GetPointer(0);
  //  pCell = static_cast<double*>(cellVolumeFraction->GetVoidPointer(0));

  // Initialize the point data to 0.
  memset(pPoint, 0, dims[0] * dims[1] * dims[2] * sizeof(double));

#ifndef NDEBUG
  // for debugging and check out of range.
  double* endPtr = pPoint + dims[0] * dims[1] * dims[2];
#endif

  //  float delProgress = (maxProgress - minProgress) / (kEnd*jEnd*iEnd) / 2;
  //  vtkIdType counter = 0;

  int index = 0;
  // Loop through the cells.
  for (k = 0; k < kEnd; ++k)
  {
    for (j = 0; j < jEnd; ++j)
    {
      for (i = 0; i < iEnd; ++i)
      {
        // if (counter % 1000 == 0 && reportProgress)
        //  {
        //  this->UpdateProgress(minProgress + delProgress*(i+j*iEnd+k*iEnd*jEnd));
        //  }
        // counter++;
        // Add cell value to all points of cell.
        double value = cellVolumeFraction->GetTuple1(index);

        assert("check: valid_range" && pPoint < endPtr);
        assert("check: valid_range" && pPoint + 1 < endPtr);
        assert("check: valid_range" && pPoint + jInc < endPtr);
        assert("check: valid_range" && pPoint + jInc + 1 < endPtr);

        *pPoint += value;
        pPoint[1] += value;
        pPoint[jInc] += value;
        pPoint[1 + jInc] += value;

        if (dimensionality == 3)
        {
          assert("check: valid_range" && pPoint + kInc < endPtr);
          assert("check: valid_range" && pPoint + kInc + 1 < endPtr);
          assert("check: valid_range" && pPoint + kInc + jInc < endPtr);
          assert("check: valid_range" && pPoint + kInc + jInc + 1 < endPtr);

          pPoint[kInc] += value;
          pPoint[kInc + 1] += value;
          pPoint[kInc + jInc] += value;
          pPoint[kInc + jInc + 1] += value;
        }

        // Increment pointers
        ++pPoint;
        ++index;
      }
      // Skip over last point to the start of the next row.
      ++pPoint;
    }
    // Skip over the last row to the start of the next plane.
    pPoint += jInc;
  }

  // Now a second pass to normalize the point values.
  // Loop through the points.
  count = 1;
  pPoint = pointVolumeFraction->GetPointer(0);

  // because we eventually modified iEnd, jEnd, kEnd to handle the
  // 2D image case, we have to recompute them.
  iEnd = dims[0] - 1;
  jEnd = dims[1] - 1;
  kEnd = dims[2] - 1;

  // counter = 0;
  for (k = 0; k <= kEnd; ++k)
  {
    // Just a fancy fast way to compute the number of cell neighbors of a
    // point.
    if (k == 1)
    {
      count = count << 1;
    }
    if (k == kEnd && kEnd > 0)
    {
      // only in 3D case, otherwise count may become zero
      // and be involved in a division by zero later on
      count = count >> 1;
    }
    for (j = 0; j <= jEnd; ++j)
    {
      // Just a fancy fast way to compute the number of cell neighbors of a
      // point.
      if (j == 1)
      {
        count = count << 1;
      }
      if (j == jEnd)
      {
        count = count >> 1;
      }
      for (i = 0; i <= iEnd; ++i)
      {
        // if (counter % 1000 == 0 && reportProgress)
        //  {
        //  this->UpdateProgress(minProgress + delProgress/2 + delProgress*(i+j*iEnd+k*iEnd*jEnd));
        //  }
        // counter++;
        // Just a fancy fast way to compute the number of cell neighbors of a
        // point.
        if (i == 1)
        {
          count = count << 1;
        }
        if (i == iEnd)
        {
          count = count >> 1;
        }
        assert("check: valid_range" && pPoint < endPtr);
        assert("check: strictly_positive_count" && count > 0);
        *pPoint = *pPoint / static_cast<double>(count);
        ++pPoint;
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkExtractCTHPart::TriggerProgressEvent(double val)
{
  double progress = this->ProgressShift + val * this->ProgressScale;
  // std::cout << "Progress: " << progress << endl;
  this->UpdateProgress(progress);
}

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

  os << indent << "VolumeArrayNames: \n";
  vtkIndent i2 = indent.GetNextIndent();
  std::vector<std::string>::iterator it;
  for (it = this->Internals->VolumeArrayNames.begin();
       it != this->Internals->VolumeArrayNames.end(); ++it)
  {
    os << i2 << *it << endl;
  }
  os << indent << "VolumeFractionSurfaceValue: " << this->VolumeFractionSurfaceValue << endl;
  os << indent << "Capping: " << this->Capping << endl;
  os << indent << "GenerateTriangles: " << this->GenerateTriangles << endl;
  os << indent << "RemoveGhostCells: " << this->RemoveGhostCells << endl;

  if (this->ClipPlane)
  {
    os << indent << "ClipPlane:\n";
    this->ClipPlane->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "ClipPlane: nullptr\n";
  }

  if (this->Controller != nullptr)
  {
    os << "Controller:" << endl;
    this->Controller->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << "No Controller." << endl;
  }
}
VTK_ABI_NAMESPACE_END