#include "vtkDiscoverRegions.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkMinimalStandardRandomSequence.h"
#include "vtkNew.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkSmartPointer.h"
#include "vtkTable.h"
#include "vtkVariant.h"
#include "vtkVector.h"
#include "vtkVectorOperators.h"
#include "vtkXMLPolyDataWriter.h"

#include "vtkObjectFactory.h"

#include "union_find.h"

#include <deque>
#include <map>
#include <set>
#include <sstream>
#include <vector>

typedef vtkVector2d vec2d;
typedef vtkVector3d vec3d;

// A edge's use of a vertex
struct VertexEdgeUse
{
  vtkIdType CellId; // index of edge in dataset
  bool Orientation; // true => face uses positive edge sense for its positive normal side. false => otherwise
  double Angle; // Angle of face (in plane with normal along edge) to first face use of edge.

  VertexEdgeUse(vtkIdType edgeId, bool orient, double angle)
    : CellId(edgeId), Orientation(orient), Angle(angle)
    {
    }
};

// Neighborhood of a vertex.
struct Vertexhood
{
  vec3d XAxis;
  vec3d BasePoint;
  std::vector<VertexEdgeUse> Bordants;
};

class Vertexhoods : public std::map<vtkIdType,Vertexhood>
{
public:
  enum
    {
    Dimension = 2
    };
};

ostream& operator << (ostream& os, Vertexhood& hood)
{
  os << " hood(" << "bpt" << hood.BasePoint <<  " x " << hood.XAxis << " Bordants [";
  for (std::vector<VertexEdgeUse>::iterator it = hood.Bordants.begin(); it != hood.Bordants.end(); ++it)
    {
    os << " cell " << it->CellId << " or " << it->Orientation << " ang " << it->Angle;
    }
  os << ")";
  return os;
};

// A face's use of an edge
struct EdgeFaceUse
{
  vtkIdType CellId; // index of face in dataset
  bool Orientation; // true => face uses positive edge sense for its positive normal side. false => otherwise
  double Angle; // Angle of face (in plane with normal along edge) to first face use of edge.
};

// Neighborhood of an edge.
struct Edgerhood
{
  vec3d Direction;
  vec3d BasePoint;
  // Axes defined by face normal of first entry in Faces.
  // They are used to keep faces ordered properly around edge:
  vec3d XAxis;
  vec3d YAxis;
  std::vector<EdgeFaceUse> Bordants;
};

typedef std::pair<vtkIdType,vtkIdType> EdgePair;

ostream& operator << (ostream& os, const EdgePair& ep)
{
  os << " edg(" << ep.first << "," << ep.second << ")";
  return os;
};

ostream& operator << (ostream& os, Edgerhood& hood)
{
  os << " hood(" << "bpt " << hood.BasePoint << " dir " << hood.Direction <<  " x " << hood.XAxis << " y " << hood.YAxis << " Bordants [";
  for (std::vector<EdgeFaceUse>::iterator it = hood.Bordants.begin(); it != hood.Bordants.end(); ++it)
    {
    os << " cell " << it->CellId << " or " << it->Orientation << " ang " << it->Angle;
    }
  os << ")";
  return os;
}

class Edgerhoods : public std::map<EdgePair,Edgerhood>
{
public:
  enum
    {
    Dimension = 3
    };
};

template<typename T>
T Add(const T& a, const T& b)
{
  T c;
  for (int i = 0; i < c.GetSize(); ++i)
    {
    c[i] = a[i] + b[i];
    }
  return c;
}

template<typename T>
T Subtract(const T& a, const T& b)
{
  T c;
  for (int i = 0; i < c.GetSize(); ++i)
    {
    c[i] = a[i] - b[i];
    }
  return c;
}

template<typename T>
void MultiplyScalar(T& a, double s)
{
  for (int i = 0; i < a.GetSize(); ++i)
    {
    a[i] *= s;
    }
}

template<typename T>
T MultiplyAdd(const T& a, const T& b, double s)
{
  T c;
  for (int i = 0; i < c.GetSize(); ++i)
    {
    c[i] = a[i] + s * b[i];
    }
  return c;
}

template<int dim>
bool EstimateNormal(vtkPoints* pts, vtkIdType npts, vtkIdType* conn, vec3d& norm);

template<>
bool EstimateNormal<2>(vtkPoints* pts, vtkIdType npts, vtkIdType* conn, vec3d& norm)
{
  if (npts < 2)
    {
    return false;
    }
  vec3d xA;
  vec3d xB;
  pts->GetPoint(conn[0], xA.GetData());
  pts->GetPoint(conn[1], xB.GetData());
  vec3d dir = xB - xA;
  vec3d z(0., 0., 1.);
  norm = dir.Cross(z);
  return norm.Normalize() < 1e-8 ? false : true;
}

template<>
bool EstimateNormal<3>(vtkPoints* pts, vtkIdType npts, vtkIdType* conn, vec3d& norm)
{
  if (npts < 3)
    {
    return false;
    }

  vec3d x[3];
  vec3d y[2];
  pts->GetPoint(conn[0], x[0].GetData());
  pts->GetPoint(conn[1], x[1].GetData());
  vtkMath::Subtract(x[1].GetData(), x[0].GetData(), y[0].GetData());
  vtkMath::Normalize(y[0].GetData());
  for (vtkIdType i = 2; i < npts; ++ i)
    {
    pts->GetPoint(conn[i], x[i % 3].GetData());
    vtkMath::Subtract(x[i % 3].GetData(), x[(i + 2) % 3].GetData(), y[(i - 1) % 2].GetData());
    vtkMath::Normalize(y[(i - 1) % 2].GetData());
    vtkMath::Cross(y[0].GetData(), y[1].GetData(), norm.GetData());
    if (vtkMath::Normalize(norm.GetData()) > 1e-8)
      {
      return true;
      }
    }
  // No sequence of 3 points along loop had a well-defined normal.
  return false;
}

static bool SpliceFaceIntoHood(
  vtkIdType faceId, vtkIdType npts, vtkIdType* conn, vtkIdType vtkNotUsed(edge), bool sense,
  Edgerhoods::iterator hood, vtkPoints* pts, double vtkNotUsed(tol))
{
  vec3d norm;
  if (!EstimateNormal<3>(pts, npts, conn, norm))
    {
    return false;
    }

  EdgeFaceUse use;
  use.CellId = faceId;
  use.Orientation = sense;
  if (!sense)
    { // flip the normal if the sense is reversed.
    MultiplyScalar(norm, -1.);
    }
  if (hood->second.Bordants.empty())
    {
    use.Angle = 0.;
    hood->second.XAxis = norm;
    vtkMath::Cross(
      hood->second.Direction.GetData(),
      hood->second.XAxis.GetData(),
      hood->second.YAxis.GetData());
    hood->second.Bordants.push_back(use);
    return true;
    }

  double x = hood->second.XAxis.Dot(norm);
  double y = hood->second.YAxis.Dot(norm);
  use.Angle = atan2(y, x);
  if (use.Angle <= 0.)
    {
    use.Angle += 2. * vtkMath::Pi();
    }
  std::vector<EdgeFaceUse>::iterator it;
  for (it = hood->second.Bordants.begin(); it != hood->second.Bordants.end(); ++it)
    {
    if (it->Angle > use.Angle)
      {
      hood->second.Bordants.insert(it, use);
      return true;
      }
    }
  // insert at end.
  hood->second.Bordants.push_back(use);
  return true;
}

int SpliceEdgeIntoVertex(
  vtkIdType edgeId, vtkIdType& idA, vtkIdType& idB, bool sense, vtkPoints* pts,
  Vertexhoods& hoods, double tol)
{
  vec3d xA;
  vec3d xB;
  pts->GetPoint(idA, xA.GetData());
  pts->GetPoint(idB, xB.GetData());
  vec3d edgeDir = xB - xA;
  std::pair<Vertexhoods::iterator,bool> srch;
  std::pair<vtkIdType,Vertexhood> entry;
  entry.first = idA;
  srch.first = hoods.find(entry.first);
  if (srch.first == hoods.end())
    {
    entry.second.BasePoint = xA;
    entry.second.XAxis = edgeDir;
    if (entry.second.XAxis.Normalize() < tol)
      { // Zero-length vector... this is a degenerate edge.
      return -1;
      }
    srch = hoods.insert(entry);
    srch.first->second.Bordants.push_back(
      VertexEdgeUse(edgeId, sense, 0.));
    }
  else
    {
    vec3d yAxis = srch.first->second.XAxis.Cross(vec3d(0,0,1));
    double angle =
      atan2(
        yAxis.Dot(edgeDir),
        srch.first->second.XAxis.Dot(edgeDir));
    VertexEdgeUse use(edgeId, sense, angle);
    std::vector<VertexEdgeUse>::iterator it;
    for (it = srch.first->second.Bordants.begin(); it != srch.first->second.Bordants.end(); ++it)
      {
      if (it->Angle > use.Angle)
        {
        srch.first->second.Bordants.insert(it, use);
        return 0;
        }
      }
    // insert at end.
    srch.first->second.Bordants.push_back(use);
    }
  return 0;
}

int SpliceCellIntoHood(
  vtkIdType edgeId, vtkIdType* conn, vtkIdType i, vtkIdType npts, vtkPoints* pts,
  Vertexhoods& hoods, double tol)
{
  vtkIdType idA = conn[i];
  vtkIdType idB = conn[(i + 1) % npts];
  // Do not insert "point" edges into vertex neighborhoods as there is
  // no way to determine their adjacent regions.
  if (idA == idB)
    {
    return -1;
    }

  // Add the edge to both of its vertex neighborhoods:
  if (SpliceEdgeIntoVertex(edgeId, idA, idB, true, pts, hoods, tol) < 0)
    {
    return -2;
    }
  if (SpliceEdgeIntoVertex(edgeId, idB, idA, false, pts, hoods, tol) < 0)
    {
    return -2;
    }
  return 0;
}

int SpliceCellIntoHood(
  vtkIdType faceId, vtkIdType* conn, vtkIdType i, vtkIdType npts, vtkPoints* pts,
  Edgerhoods& hoods, double tol)
{
  vtkIdType idA = conn[i];
  vtkIdType idB = conn[(i + 1) % npts];
  // Ignore topological degeneracies if we have enough vertices to continue:
  if (idA == idB)
    {
    return -1;
    }

  // Look up the edge by its canonical order and remember
  // its orientation (sense) relative to that order.
  std::pair<EdgePair,Edgerhood> entry;
  bool sense = idA < idB; // true -> edge and face are oriented consistently
  entry.first.first = sense ? idA : idB;
  entry.first.second = sense ? idB : idA;
  std::pair<Edgerhoods::iterator,bool> srch;
  srch.first = hoods.find(entry.first);
  // If we haven't seen the edge before, initialize a hood for it:
  if (srch.first == hoods.end())
    {
    srch = hoods.insert(entry);
    vec3d xA;
    vec3d xB;
    pts->GetPoint(entry.first.first, xA.GetData());
    pts->GetPoint(entry.first.second, xB.GetData());
    srch.first->second.BasePoint = xA;
    srch.first->second.Direction = Subtract(xB, xA);
    // If the edge is too short (i.e., pts[idA] == pts[idB]),
    // treat it as a single point as long as we have enough
    // other vertices. Of course, if the edge is also non-manifold
    // or bordant to 2 regions, we could be in trouble.
    if (srch.first->second.Direction.Normalize() < tol)
      {
      hoods.erase(srch.first);
      return -1;
      }
    }
  if (!SpliceFaceIntoHood(faceId, npts, conn, i, sense, srch.first, pts, tol))
    { // Error out as splicing only fails when the whole face is degenerate.
    return -2;
    }
  return 0;
}

template<typename H>
static void BuildNeighborhoods(vtkPolyData* pd, H& hoods)
{
  int dim = H::Dimension;
  hoods.clear();
  vtkCellArray* cells = (dim == 2 ? pd->GetLines() : pd->GetPolys());
  vtkPoints* pts = pd->GetPoints();
  if (!cells || !pts)
    {
    return;
    }

  // Compute tolerance for when distances are considered degenerate:
  double diam = pd->GetLength();
  double tol = diam * VTK_DBL_EPSILON * 2.;
  tol = (tol <= 0. ? VTK_DBL_EPSILON : tol);

  // Now traverse all (dim-1)-boundaries of all cells and insert the cell into
  // each (dim-1)-neighborhood (vertex neighborhoods for edge cells in 2-D,
  // edge neighborhoods for face cells in 3-D).
  vtkIdType npts;
  vtkIdType* conn;
  vtkIdType cellId =
    (dim == 2 ?
     pd->GetNumberOfVerts() :
     pd->GetNumberOfVerts() + pd->GetNumberOfLines());
  for (cells->InitTraversal(); cells->GetNextCell(npts, conn); ++cellId)
    {
    vtkIdType distinctVertsCriterion = npts;
    if (npts < dim)
      continue;
    // Iterate over the (dim-1)-boundaries of the cell.
    // Polygons (dim 3) are implicitly closed loops (the last edge
    // connecting the start and end vertices is implied).
    // Polylines (dim 2) may be open or closed, but are always explicit.
    vtkIdType loopBounds = (dim == 2 ? npts - 1 : npts);
    for (vtkIdType i = 0; i < loopBounds; ++i)
      {
      if (distinctVertsCriterion < dim)
        {
        break;
        }
      int status = SpliceCellIntoHood(cellId, conn, i, npts, pts, hoods, tol);
      switch (status)
        {
      case 0:
        // OK. Do nothing.
        break;
      case -1:
        // Degenerate coordinates. Skip.
        --distinctVertsCriterion;
        break;
      case -2:
        // Degenerate cell. Stop.
        distinctVertsCriterion = -1;
        break;
        }
      }
    if (distinctVertsCriterion < dim)
      {
      vtkGenericWarningMacro(
        << "Cell " << cellId << " has only "
        << distinctVertsCriterion << " distinct vertices");
      }
    }
}

namespace {

template<typename N>
struct RegionTracker
{
  // Each entry in Sets represents the region to one side of one triangle (or model face).
  UnionFind Sets;
  // Region to each side of each triangle (or model face when ModelFaceArrayName is set).
  vtkSmartPointer<vtkIdTypeArray> ModelRegions;
  // Region to each side of each *proper* model face (bounding 2 or less regions)
  vtkSmartPointer<vtkIdTypeArray> ReconciledModelRegions;
  // Map from model-facet ID to a tuple in ModelRegions
  std::map<vtkIdType,vtkIdType> ModelMap;
  // Renumbering of Sets IDs to sequential integers after merge along shared edges is complete.
  std::map<vtkIdType,vtkIdType> Collapse;
  // Faces in counterclockwise order around each edge.
  N CellNeighborhoods;
  // ID for each region assigned by points on input port 3.
  vtkNew<vtkIntArray> RegionIds;
  // ID of region that is exterior to all others
  vtkIdType ExteriorId;
  // Arrays to be associated with each identified region.
  vtkNew<vtkPointData> RegionAttributes;
  // A table describing containment relationships between connected components.
  vtkNew<vtkTable> ContainmentRelationships;
  // A column in ContainmentRelationships whose values are region IDs of contained shells
  vtkNew<vtkIdTypeArray> ContainedShellIds;
  // A column in ContainmentRelationships whose values are cells on contained shells
  vtkNew<vtkIdTypeArray> ContainedShellCells;
  // A column in ContainmentRelationships whose values are cells on contained shells
  vtkNew<vtkIntArray> ContainedShellSense;
  // A column in ContainmentRelationships whose values are region IDs of containers (or -1)
  vtkNew<vtkIdTypeArray> ContainerShellIds;
  // A column in ContainmentRelationships whose values are cells on containing shells (or -1)
  vtkNew<vtkIdTypeArray> ContainerShellCells;
  // A column in ContainmentRelationships whose values are cells on containing shells (or -1)
  vtkNew<vtkIntArray> ContainerShellSense;
  // An optional polydata pointer in which a point should be stored for each region
  vtkPolyData* RegionPoints;

  RegionTracker()
    {
    this->ExteriorId = -1;
    this->ContainedShellIds->SetName("ContainedShellIds");
    this->ContainedShellCells->SetName("ContainedShellCells");
    this->ContainedShellSense->SetName("ContainedShellSense");
    this->ContainerShellIds->SetName("ContainerShellIds");
    this->ContainerShellCells->SetName("ContainerShellCells");
    this->ContainerShellSense->SetName("ContainerShellSense");
    this->ContainmentRelationships->AddColumn(this->ContainedShellIds.GetPointer());
    this->ContainmentRelationships->AddColumn(this->ContainedShellCells.GetPointer());
    this->ContainmentRelationships->AddColumn(this->ContainedShellSense.GetPointer());
    this->ContainmentRelationships->AddColumn(this->ContainerShellIds.GetPointer());
    this->ContainmentRelationships->AddColumn(this->ContainerShellCells.GetPointer());
    this->ContainmentRelationships->AddColumn(this->ContainerShellSense.GetPointer());
    this->RegionPoints = NULL;
    }

  // Keep a record of which connected components are contained in which loops.
  void AddContainmentRelation(
    vtkIdType cellOnShell, bool cellSense,
    vtkIdType cellOnContainer, bool containerSense)
    {
    this->ContainedShellIds->InsertNextValue(-1);
    this->ContainedShellCells->InsertNextValue(cellOnShell);
    this->ContainedShellSense->InsertNextValue(cellSense);
    this->ContainerShellIds->InsertNextValue(-1);
    this->ContainerShellCells->InsertNextValue(cellOnContainer);
    this->ContainerShellSense->InsertNextValue(containerSense);
    }

  // Map region IDs in ContainmentRelationships table to their collapsed IDs
  void CollapseContainmentRelationships()
    {
    vtkIdType numRelations = this->ContainedShellIds->GetNumberOfTuples();
    vtkIdType* containedRegions = this->ContainedShellIds->GetPointer(0);
    vtkIdType* containerRegions = this->ContainerShellIds->GetPointer(0);
    vtkIdType* containedCellIds = this->ContainedShellCells->GetPointer(0);
    vtkIdType* containerCellIds = this->ContainerShellCells->GetPointer(0);
    int* containedSense = this->ContainedShellSense->GetPointer(0);
    int* containerSense = this->ContainerShellSense->GetPointer(0);
    for (vtkIdType i = 0; i < numRelations; ++i)
      {
      /*
      cout
        << "Contained " << containedCellIds[i] << "," << containedSense[i]
        << "," << this->Collapse[this->Sets.Find(this->ModelRegions->GetValue(3 * containedCellIds[i] + (containedSense[i] ? 2 : 1)))]
        << " container " << containerCellIds[i] << "," << containerSense[i]
        << "," << this->Collapse[this->Sets.Find(this->ModelRegions->GetValue(3 * containerCellIds[i] + (containerSense[i] ? 2 : 1)))]
        << "\n";
        */
      containedRegions[i] = this->Collapse[
        this->Sets.Find(
          this->ModelRegions->GetValue(
            3 * containedCellIds[i] + (containedSense[i] ? 2 : 1)))];
      containerRegions[i] = this->Collapse[
        this->Sets.Find(
          this->ModelRegions->GetValue(
            3 * containerCellIds[i] + (containerSense[i] ? 2 : 1)))];
      }
    }

  void DumpHoods()
    {
    typename N::iterator it;
    for (it = this->CellNeighborhoods.begin(); it != this->CellNeighborhoods.end(); ++it)
      {
      cout << it->first << ": " << it->second << "\n";
      }
    }

  void Dump()
    {
    vtkIdType entry[3];
    for (
      std::map<vtkIdType,vtkIdType>::iterator it = this->ModelMap.begin();
      it != this->ModelMap.end();
      ++it)
      {
      this->ModelRegions->GetTupleValue(it->second, entry);
      vtkIdType r0 = this->Sets.Find(entry[1]);
      vtkIdType r1 = this->Sets.Find(entry[2]);
      cout
        << "  Face " << it->first << " facet " << it->second
        << " bounds regions " << r0 << "," << r1 << "\n";
      }
    }

  void DumpCollapsed()
    {
    vtkIdType entry[3];
    for (
      std::map<vtkIdType,vtkIdType>::iterator it = this->ModelMap.begin();
      it != this->ModelMap.end();
      ++it)
      {
      this->ModelRegions->GetTupleValue(it->second, entry);
      vtkIdType r0 = this->Collapse[this->Sets.Find(entry[1])];
      vtkIdType r1 = this->Collapse[this->Sets.Find(entry[2])];
      cout
        << "  Face " << it->first << " facet " << it->second
        << " bounds regions " << r0 << "," << r1 << "\n";
      }
    }
};

} // namespace

// -----------------------------------------------------------------------------
vtkStandardNewMacro(vtkDiscoverRegions);

vtkDiscoverRegions::vtkDiscoverRegions()
{
  this->FaceGroupArrayName = NULL;
  this->ModelFaceArrayName = NULL;
  this->RegionGroupArrayName = NULL;
  this->SetFaceGroupArrayName("FaceGroups");
  this->ReportRegionsByModelFace = 0;
  this->GenerateRegionInteriorPoints = 0;
  // port 0: triangulated facets
  // port 1: volume hole points
  // port 2: region assignment points
  this->SetNumberOfInputPorts(3);
  // port 0: input marked with region information
  // port 1: table of containment relations between shells
  // port 2: points for discovered regions
  this->SetNumberOfOutputPorts(3);
}

vtkDiscoverRegions::~vtkDiscoverRegions()
{
  this->SetFaceGroupArrayName(NULL);
  this->SetModelFaceArrayName(NULL);
  this->SetRegionGroupArrayName(NULL);
}

void vtkDiscoverRegions::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "ReportRegionsByModelFace: " << this->ReportRegionsByModelFace << "\n";
  os << indent << "GenerateRegionInteriorPoints: " << this->GenerateRegionInteriorPoints << "\n";
  os << indent << "FaceGroupArrayName: " << this->FaceGroupArrayName << "\n";
  os << indent << "ModelFaceArrayName: " << this->ModelFaceArrayName << "\n";
  os << indent << "RegionGroupArrayName: " << this->RegionGroupArrayName << "\n";
}

template<typename T>
struct vtkDiscoverRegionsFacetFromCellArray
{
  vtkIdType operator () (vtkIdType poly) const
    {
    return this->Array->GetValue(poly);
    }
  std::set<vtkIdType> FacetIds(vtkPolyData*) const
    {
    std::set<vtkIdType> uniques;
    for (vtkIdType i = 0; i <= this->Array->GetMaxId(); ++i)
      {
      uniques.insert(this->Array->GetValue(i));
      }
    return uniques;
    }
  vtkIdType CellForFacet(vtkIdType facet)
    {
    for (vtkIdType i = 0; i <= this->Array->GetMaxId(); ++i)
      {
      if (facet == this->Array->GetValue(i))
        {
        return i;
        }
      }
    return -1;
    }
  vtkIdTypeArray* GetFaceGroupArray()
    {
    return this->Array;
    }
  bool AreOversubscribedFacetsPossible(vtkIdTypeArray* other) const
    {
    return other != this->Array;
    }
  T* Array;
};

struct vtkDiscoverRegionsOneFacetPerCell
{
  vtkIdType operator () (vtkIdType poly) const
    {
    return poly;
    }
  std::set<vtkIdType> FacetIds(vtkPolyData* src) const
    {
    std::set<vtkIdType> uniques;
    for (vtkIdType i = 0; i < src->GetNumberOfCells(); ++i)
      {
      uniques.insert(i);
      }
    return uniques;
    }
  vtkIdType CellForFacet(vtkIdType facet)
    {
    return facet; // facet == cell
    }
  vtkIdTypeArray* GetFaceGroupArray()
    {
    return NULL;
    }
  bool AreOversubscribedFacetsPossible(vtkIdTypeArray*) const
    {
    return true;
    }
};

template<typename T, typename N>
void InitializeRegions(vtkPolyData* surface, RegionTracker<N>& regions, const T& cell2facet)
{
  std::set<vtkIdType> modelFacets = cell2facet.FacetIds(surface);
  regions.ModelRegions = vtkSmartPointer<vtkIdTypeArray>::New();
  regions.ModelRegions->SetName("ModelFaceRegionsMap");
  regions.ModelRegions->SetNumberOfComponents(3);
  regions.ModelRegions->SetComponentName(0, "ModelFace");
  regions.ModelRegions->SetComponentName(1, "BackfaceRegion");
  regions.ModelRegions->SetComponentName(2, "FrontfaceRegion");
  regions.ModelRegions->SetNumberOfTuples(modelFacets.size());
  vtkIdType i = 0;
  for (std::set<vtkIdType>::iterator it = modelFacets.begin(); it != modelFacets.end(); ++it, ++i)
    {
    vtkIdType tuple[3] = { *it, regions.Sets.NewSet(), regions.Sets.NewSet() };
    regions.ModelRegions->SetTupleValue(i, tuple);
    regions.ModelMap[*it] = i;
    }
}

template<typename T, typename N>
void MergeRegions(vtkPolyData* pdIn, RegionTracker<N>& regions, const T& cell2facet)
{
  pdIn->BuildCells();
  BuildNeighborhoods<N>(pdIn, regions.CellNeighborhoods);
  // Iterate over each edge, E = {e0,e1} with attached faces {f_i = 0,n}
  typename N::iterator it;
  for (
    it = regions.CellNeighborhoods.begin();
    it != regions.CellNeighborhoods.end();
    ++it)
    {
    //cout << "Considering neighborhood of " << it->first << "\n";
    // Faces are orderedCCW along positive edge dir.
    vtkIdType numFaces = static_cast<vtkIdType>(it->second.Bordants.size());
    for (vtkIdType f = 0; f < numFaces; ++f)
      {
      vtkIdType cellIdA = it->second.Bordants[f].CellId;
      vtkIdType cellIdB = it->second.Bordants[(f + 1) % numFaces].CellId;

      vtkIdType modelRegionA = regions.ModelMap[cell2facet(cellIdA)];
      vtkIdType modelRegionB = regions.ModelMap[cell2facet(cellIdB)];

      vtkIdType materialA =
        regions.ModelRegions->GetValue(modelRegionA * 3 +
          (it->second.Bordants[f].Orientation ? 2 : 1));
      vtkIdType materialB =
        regions.ModelRegions->GetValue(modelRegionB * 3 +
          (it->second.Bordants[(f + 1) % numFaces].Orientation ? 1 : 2));
      /*
      cout << "  merge "
        << "cell " << cellIdA << " reg " << materialA << " with "
        << "cell " << cellIdB << " reg " << materialB << "\n";
        */

      regions.Sets.MergeSets(materialA, materialB);
      }
    }
  //regions.Dump();
}

struct RayHitRecord
{
  vec3d S;
  double T;
  vec3d X;
  vtkIdType RegionInfo[3];
  vtkIdType CellId;
  bool Sense;

  bool operator < (const RayHitRecord& other) const
    { return this->T < other.T; }
};

typedef std::multiset<RayHitRecord> RayHitRecords;

template<typename N>
bool NearCorner(const vtkVector3d& param)
{
  int dim = N::Dimension;
  bool interior = true;
  for (int i = 0; i < dim; ++i)
    {
    // We are "well into" the interior if not near param[i] == 0 or param[i] == 1:
    interior &= (fabs(param[i]) > 1e-8 && fabs(1. - param[i]) > 1e-8);
    }
  return !interior;
};

template<typename T, typename N>
bool IntersectWithRay(
  vtkPolyData* pdIn, RegionTracker<N>& regions, T& cell2facet,
  vec3d& basePoint, const vec3d& direction,
  RayHitRecords& hits, vtkIdType except = -1)
{
  bool allHitsOk = true;
  int dim = N::Dimension;
  hits.clear();
  RayHitRecord hit;
  vtkPoints* pts = pdIn->GetPoints();
  double diam = pdIn->GetLength();
  vec3d p2 = MultiplyAdd(basePoint, direction, diam);
  vtkIdType begin = pdIn->GetNumberOfVerts() + (dim == 2 ? 0 : pdIn->GetNumberOfLines());
  vtkIdType end = begin + (dim == 2 ? pdIn->GetNumberOfLines() : pdIn->GetNumberOfPolys());
  for (hit.CellId = begin; hit.CellId < end; ++hit.CellId)
    {
    if (hit.CellId == except)
      { // Don't intersect with a face we're told to ignore.
      continue;
      }
    vtkCell* cell = pdIn->GetCell(hit.CellId);
    int dummySubId;
    if (
      cell->IntersectWithLine(
        basePoint.GetData(), p2.GetData(), 1e-8 * diam,
        hit.T, hit.X.GetData(), hit.S.GetData(), dummySubId))
      {
      vtkIdType npts;
      vtkIdType* conn;
      vec3d norm;
      pdIn->GetCellPoints(hit.CellId, npts, conn);
      EstimateNormal<N::Dimension>(pts, npts, conn, norm);
      vtkIdType modelFacet = regions.ModelMap[cell2facet(hit.CellId)];
      regions.ModelRegions->GetTupleValue(modelFacet, hit.RegionInfo);
      double normDotDir;
      if ((normDotDir = norm.Dot(direction)) < 0.)
        {
        vtkIdType tmp = hit.RegionInfo[1];
        hit.RegionInfo[1] = hit.RegionInfo[2];
        hit.RegionInfo[2] = tmp;
        hit.Sense = false;
        }
      else
        {
        hit.Sense = true;
        }
      allHitsOk |= !(fabs(normDotDir) < 1e-8 || NearCorner<N>(hit.S));
      for (int j = 1; j < 3; ++j)
        {
        hit.RegionInfo[j] = regions.Sets.Find(hit.RegionInfo[j]);
        }
      hits.insert(hit);
      }
    }
  return true; //allHitsOk;
}

struct RegionHitCountRecord
{
  double FirstHit;
  bool DoesFirstHitExit; // is the first hit a region entry(F) or exit(T) event?
  int Count; // sum of signed entry(+1) + exit(-1) events

  RegionHitCountRecord()
    {
    this->FirstHit = vtkMath::Inf();
    this->DoesFirstHitExit = false;
    this->Count = 0;
    }
};

static bool ContainingShellFromHits(
  RayHitRecords& hits, vtkIdType& region, vtkIdType& hitCell, bool& hitSense, bool returnInwardShell)
{
  //cout << "--\n";
  std::deque<RayHitRecords::iterator> inwardHits;
  std::deque<RayHitRecords::iterator> outwardHits;
  for (RayHitRecords::iterator it = hits.begin(); it != hits.end(); ++it)
    {
    /*
    cout
      << "    * " << it->T << "  x" << it->X << " modfacet " << it->RegionInfo[0]
      << " regs: - " << it->RegionInfo[1] << " + " << it->RegionInfo[2] << "\n";
      */

    if (!outwardHits.empty() && outwardHits.front()->RegionInfo[2] == it->RegionInfo[1])
      {
      outwardHits.pop_front();
      }
    else
      {
      inwardHits.push_front(it);
      }
    if (!inwardHits.empty() && inwardHits.front()->RegionInfo[1] == it->RegionInfo[2])
      {
      inwardHits.pop_front();
      }
    else
      {
      outwardHits.push_front(it);
      }
    }
  if (returnInwardShell)
    {
    if (inwardHits.empty())
      { // Didn't hit anything but air.
      region = -1;
      hitCell = -1;
      //cout << "  REGION -1\n";
      return false;
      }
    region = inwardHits.back()->RegionInfo[1];
    hitCell = inwardHits.back()->CellId;
    hitSense = ! inwardHits.back()->Sense;
    //cout << "  REGION " << inwardHits.back()->RegionInfo[1] << "\n";
    }
  else
    {
    //if (outwardStack.empty())
    if (outwardHits.empty())
      { // Didn't hit anything but air.
      region = -1;
      hitCell = -1;
      //cout << "  REGION -1\n";
      return false;
      }
    region = outwardHits.back()->RegionInfo[2];
    hitCell = outwardHits.back()->CellId;
    hitSense = outwardHits.back()->Sense;
    //cout << "  REGION " << outwardHits.back()->RegionInfo[2] << "\n";
    }
  return true;
}

static bool ContainingShellFromHits(
  RayHitRecords& hits, vtkIdType& region, bool returnInwardShell = true)
{
  vtkIdType dummyCell;
  bool dummySense;
  return ContainingShellFromHits(hits, region, dummyCell, dummySense, returnInwardShell);
}

template<typename T, typename N>
void FindPointsInRegions(
  vtkPolyData* pdIn, RegionTracker<N>& regions, T& cell2facet,
  vtkMinimalStandardRandomSequence* rnd)
{
  std::set<vtkIdType> shells = regions.Sets.Roots();
  std::set<vtkIdType>::iterator shellIt;
  std::map<vtkIdType,vtkVector3d> pointsByRegion;

  vtkPoints* pts = pdIn->GetPoints();
  vtkIdType npts; // number of points in one cell (not pts)
  vtkIdType* conn; // cell connectivity

  // Iterate over the shells
  for (
    shellIt = shells.begin();
    shellIt != shells.end();
    ++shellIt)
    {
    if (pointsByRegion.find(*shellIt) != pointsByRegion.end())
      {
      continue; // Another ray already provided a point in the region.
      }
    //cout << "\nPoints for Shell " << *shellIt << "\n";

    // Note that shell IDs have a special relationship to ModelMap;
    // they are assigned to ModelMap in sequential, 0-based pairs.
    // So, given a shell ID, we can divide by two and know an offset
    // into ModelMap. The remainder indicates whether the negative
    // (1 remainder) or positive (0 remainder) orientation of facets
    // correspond to the shell ID. From this offset into ModelMap,
    // we can then find the Facet ID for a shell and from there obtain
    // a cell on the facet for the shell.
    bool sense = (*shellIt) % 2 ? true : false; // positive orientation?
    vtkIdType regionInfo[3];
    regions.ModelRegions->GetTupleValue(*shellIt / 2, regionInfo);
    vtkIdType cellOnShell = cell2facet.CellForFacet(regionInfo[0]);

    // Get cell normal
    vec3d norm;
    pdIn->GetCellPoints(cellOnShell, npts, conn);
    EstimateNormal<N::Dimension>(pts, npts, conn, norm);

    // Flip normal if shell and cell sense are opposite
    if (!sense)
      {
      MultiplyScalar(norm, -1.);
      }

    // Now we have a cell on the shell and know the shell normal.
    // Shoot rays "near" the normal to find which shells contain us.
    // TODO:
    //  1. Ideally here we would loop over several different pcenter
    //     coordinates for numerical stability
    //  2. This assumes we are dealing with a cell whose average
    //     coordinate is in its interior
    vec3d basePt(0.);
    vec3d tmp;
    double invNpts = 1. / npts;
    for (int i = 0; i < npts; ++i)
      {
      //cout << " " << conn[i];
      pts->GetPoint(conn[i], tmp.GetData());
      basePt = MultiplyAdd(basePt, tmp, invNpts);
      }
    bool rayHitsAllOk = false; // Did the ray hit any cell vertices or run tangent near a cell?
    std::map<vtkIdType, double> containerCount;
    int numTries = 0;
    // Perturb the normal at least 5 times and repeat until
    // over half the container answers agree. Or we've tried 15
    // times, which is a local approximation to an infinite loop.
    while (!rayHitsAllOk && numTries < 15)
      {
      // Perturb the unit normal by at most a half-unit vector.
      // This prevents the normal from disappearing, or
      // pointing into the solid, or along a tangent to the solid.
      vec3d pnorm(0.);
      for (int i = 0; i < N::Dimension; ++i, rnd->Next())
        {
        pnorm[i] = rnd->GetValue();
        }
      double scale = pnorm.SquaredNorm();
      // Renormalize after adding.
      (pnorm = MultiplyAdd(norm, pnorm, scale > 0.25 ? 0.5 / sqrt(scale) : 1.)).Normalize();
      /*
      cout
        << "  Tries " << numTries << " Cell " << cellOnShell << " np " << npts << " ["
        << "] sense " << (sense ? "+" : "-") << " pt " << basePt << " dir " << pnorm
        << "\n";
        */
      RayHitRecords hits;
      rayHitsAllOk |= IntersectWithRay(pdIn, regions, cell2facet, basePt, pnorm, hits, -1);
      if (rayHitsAllOk)
        {
        RayHitRecords::iterator i0 = hits.begin();
        RayHitRecords::iterator i1 = i0;
        rayHitsAllOk = false;
        if (i0 != hits.end())
          {
          for (++i1; i1 != hits.end(); ++i0, ++i1)
            {
            vtkVector3d pt = 0.5 * (i0->X + i1->X);
            int region = FindRegionContainingPoint(pdIn, pt, regions, cell2facet, rnd);
            //cout << "  Region " << region << "  " << pt << "\n";
            if (region == *shellIt)
              { // We got the region we need
              rayHitsAllOk = true;
              }
            if (pointsByRegion.find(region) == pointsByRegion.end() && region == *shellIt)
              {
              pointsByRegion[region] = pt;
              }
            }
          }
        }
      ++numTries;
      }
    }
  vtkNew<vtkPoints> rpts;
  vtkNew<vtkIdTypeArray> regionId;
  regionId->SetName("Region");
  regionId->SetNumberOfTuples(pointsByRegion.size());
  rpts->SetNumberOfPoints(pointsByRegion.size());
  regions.RegionPoints->SetPoints(rpts.GetPointer());
  regions.RegionPoints->GetPointData()->SetScalars(regionId.GetPointer());
  //cout << "Points in regions: (" << shells.size() << " == " << pointsByRegion.size() << "?)\n";
  std::map<vtkIdType,vtkVector3d>::iterator pit;
  vtkIdType i = 0;
  for (pit = pointsByRegion.begin(); pit != pointsByRegion.end(); ++pit, ++i)
    {
    rpts->SetPoint(i, pit->second.GetData());
    regionId->SetValue(i, pit->first);
    //cout << "  " << pit->first << "   " << pit->second << "\n";
    }
}

template<typename T, typename N>
vtkIdType FindRegionContainingPoint(
  vtkPolyData* pdIn, vtkVector3d& pt,
  RegionTracker<N>& regions, T& cell2facet,
  vtkMinimalStandardRandomSequence* rnd)
{
  // Choose a random direction:
  vec3d r;
  do
    {
    for (int i = 0; i < N::Dimension; ++i, rnd->Next())
      {
      r[i] = rnd->GetValue();
      }
    }
  while (r.Normalize() < 1e-5);
  // Find shell intersections
  RayHitRecords hits;
  IntersectWithRay(pdIn, regions, cell2facet, pt, r, hits);
  vtkIdType region;
  if (!ContainingShellFromHits(hits, region) || region < 0)
    {
    region = regions.ExteriorId;
    }
  return region;
}

template<typename T, typename N>
vtkIdType DiscoverNestings(
  vtkPolyData* pdIn, RegionTracker<N>& regions, T& cell2facet)
{
  // For each *unique* entry in regions.Sets,
  //   Choose a random point, P, on a seed face, F
  //   For each co-face of F,
  //     Intersect a ray (from P normal to the co-face) with other faces not in the same region
  //     Choose the smallest t > 0 intersection (say with face I)
  //     If the intersection is close to an edge or vertex of I, choose a different face F and retry
  //     Resolve the region ID from the intersected face I (based on the orientation of I relative to P)
  //     Merge the regions for the co-face of F and the co-face of I that face each other.

  vtkPoints* pts = pdIn->GetPoints();
  vtkIdType npts; // number of points in one cell (not pts)
  vtkIdType* conn; // cell connectivity

  // Get a list of shells:
  std::set<vtkIdType> shells = regions.Sets.Roots();
  std::set<vtkIdType>::iterator shellIt;

  // A source of random vectors for repeated trials
  vtkNew<vtkMinimalStandardRandomSequence> rnd;

  // If requested by user, compute a point in the interior of each region:
  if (regions.RegionPoints)
    {
    FindPointsInRegions(pdIn, regions, cell2facet, rnd.GetPointer());
    }

  // Iterate over the shells
  for (
    shellIt = shells.begin();
    shellIt != shells.end();
    ++shellIt)
    {
    //cout << "\nNesting for Shell " << *shellIt << "\n";
    // Note that shell IDs have a special relationship to ModelMap;
    // they are assigned to ModelMap in sequential, 0-based pairs.
    // So, given a shell ID, we can divide by two and know an offset
    // into ModelMap. The remainder indicates whether the negative
    // (1 remainder) or positive (0 remainder) orientation of facets
    // correspond to the shell ID. From this offset into ModelMap,
    // we can then find the Facet ID for a shell and from there obtain
    // a cell on the facet for the shell.
    bool sense = (*shellIt) % 2 ? true : false; // positive orientation?
    vtkIdType regionInfo[3];
    regions.ModelRegions->GetTupleValue(*shellIt / 2, regionInfo);
    vtkIdType cellOnShell = cell2facet.CellForFacet(regionInfo[0]);

    // Get cell normal
    vec3d norm;
    pdIn->GetCellPoints(cellOnShell, npts, conn);
    EstimateNormal<N::Dimension>(pts, npts, conn, norm);

    // Flip normal if shell and cell sense are opposite
    if (!sense)
      {
      MultiplyScalar(norm, -1.);
      }

    // Now we have a cell on the shell and know the shell normal.
    // Shoot rays "near" the normal to find which shells contain us.
    // TODO:
    //  1. Ideally here we would loop over several different pcenter
    //     coordinates for numerical stability
    //  2. This assumes we are dealing with a cell whose average
    //     coordinate is in its interior
    vec3d basePt(0.);
    vec3d tmp;
    double invNpts = 1. / npts;
    for (int i = 0; i < npts; ++i)
      {
      //cout << " " << conn[i];
      pts->GetPoint(conn[i], tmp.GetData());
      basePt = MultiplyAdd(basePt, tmp, invNpts);
      }
    std::map<vtkIdType, double> containerCount;
    vtkIdType container = -1; // The region ID of the innermost container hit by a test ray.
    vtkIdType containerCell = -1; // The ID of the cell hit on the container by a test ray.
    bool containerSense = false; // Whether the containerCell was hit on its positive (true) or negative (false) side.
    int numTries = 0;
    // Perturb the normal at least 5 times and repeat until
    // over half the container answers agree. Or we've tried 15
    // times, which is a local approximation to an infinite loop.
    while (
      numTries < 5 || (
        (!containerCount.empty() && containerCount[container] / numTries < 0.5) &&
        numTries < 15))
      {
      // Perturb the unit normal by at most a half-unit vector.
      // This prevents the normal from disappearing, or
      // pointing into the solid, or along a tangent to the solid.
      vec3d pnorm(0.);
      for (int i = 0; i < N::Dimension; ++i, rnd->Next())
        {
        pnorm[i] = rnd->GetValue();
        }
      double scale = pnorm.SquaredNorm();
      // Renormalize after adding.
      (pnorm = MultiplyAdd(norm, pnorm, scale > 0.25 ? 0.5 / sqrt(scale) : 1.)).Normalize();
      /*
      cout
        << "  Tries " << numTries << " Cell " << cellOnShell << " np " << npts << " ["
        << "] sense " << (sense ? "+" : "-") << " pt " << basePt << " dir " << pnorm
        << " container " << container
        << " ok " << (containerCount.empty() ? -1. : containerCount[container] / numTries)
        << "\n";
        */
      RayHitRecords hits;
      IntersectWithRay(pdIn, regions, cell2facet, basePt, pnorm, hits, cellOnShell);
      vtkIdType tmpCon;
      vtkIdType tmpCell;
      bool tmpSense;
      if (ContainingShellFromHits(hits, tmpCon, tmpCell, tmpSense, true))
        {
        containerCount[tmpCon] += 1.;
        if (containerCount[tmpCon] > containerCount[container])
          {
          container = tmpCon;
          containerCell = tmpCell;
          containerSense = tmpSense;
          }
        }

      ++numTries;
      }
    /*
    cout
      << "  Tries " << numTries << " Cell " << cellOnShell << " np " << npts << " ["
      << "] sense " << (sense ? "+" : "-") << " pt " << basePt << " dir " << norm
      << " container " << container << " ok " << (containerCount.empty() ? "F" : "T")
      << "\n";
      */
    vtkIdType actualShellId = regions.Sets.Find(*shellIt);
    if (!containerCount.empty() && container >= 0)
      {
      if (actualShellId == container)
        {
        //cout << "Shell " << container << " stabbed itself\n";
        }
      else
        {
        /*
        cout
          << "Shell " << *shellIt << " cell " << cellOnShell << " sense " << sense
          << " inside " << container << " cell " << containerCell << " sense " << containerSense << "\n";
          */
        regions.AddContainmentRelation(cellOnShell, sense, containerCell, containerSense);
        regions.Sets.MergeSets(*shellIt, container);
        }
      }
    else
      {
      //cout << "Shell " << *shellIt << " is exterior\n";
      if (regions.ExteriorId < 0)
        {
        regions.ExteriorId = *shellIt;
        }
      else
        {
        // This is not really a containment relationship, so do not
        // call regions.AddContainmentRelation(...); if you change
        // this, then the cell and sense for the exterior region
        // must be valid.
        regions.Sets.MergeSets(*shellIt, regions.ExteriorId);
        }
      }
    }
  return regions.ExteriorId;
}

struct ModelFacetSplitRec
{
  vtkIdType OriginalFacet;
  vtkIdType Regions[2];

  ModelFacetSplitRec(vtkIdType orig, vtkIdType r0, vtkIdType r1)
    {
    this->OriginalFacet = orig;
    this->Regions[0] = r0 < r1 ? r0 : r1;
    this->Regions[1] = r0 < r1 ? r1 : r0;
    }

  ModelFacetSplitRec(const ModelFacetSplitRec& other)
    {
    this->OriginalFacet = other.OriginalFacet;
    for (int i = 0; i < 2; ++i)
      {
      this->Regions[i] = other.Regions[i];
      }
    }

  bool operator < (const ModelFacetSplitRec& other) const
    {
    return
      (this->OriginalFacet < other.OriginalFacet ? true :
       (this->OriginalFacet > other.OriginalFacet ? false :
        (this->Regions[0] < other.Regions[0] ? true :
         (this->Regions[0] > other.Regions[0] ? false :
          (this->Regions[1] < other.Regions[1] ? true :
           (this->Regions[1] > other.Regions[1] ? false :
            false))))));
    }
};

class ModelFacetSplitMap : public std::map<ModelFacetSplitRec,vtkIdType>
{
public:
  ModelFacetSplitMap(vtkIdType numOriginalModelFacets)
    {
    this->NextModelFacet = numOriginalModelFacets;
    }
  vtkIdType SplitFaceByRegion(vtkIdType orig, vtkIdType r0, vtkIdType r1)
    {
    std::pair<ModelFacetSplitRec,vtkIdType>
      rec(ModelFacetSplitRec(orig, r0, r1), this->NextModelFacet);
    iterator it = this->find(rec.first);
    if (it == this->end())
      { // Insert a new model facet
      if (this->Originals.find(orig) == this->Originals.end())
        {
        rec.second = orig;
        this->Originals.insert(orig);
        }
      else
        {
        ++this->NextModelFacet;
        }
      it = this->insert(rec).first;
      }
    return it->second;
    }
  vtkIdType NextModelFacet;
  std::set<vtkIdType> Originals;
};

template<typename T, typename N>
void ReconcileModelFacets(
  vtkPolyData* pdIn, const char* modelFaceArrayName,
  RegionTracker<N>& regions, T& cell2facet)
{
  vtkIdTypeArray* modelFaceArray =
    vtkIdTypeArray::SafeDownCast(
      pdIn->GetCellData()->GetArray(modelFaceArrayName));
  if (!modelFaceArray)
    {
    vtkGenericWarningMacro("No array \"" << modelFaceArrayName << "\" found.");
    return;
    }
  vtkIdType numModelFaces = modelFaceArray->GetValueRange(0)[1];
  if (modelFaceArray->GetValueRange(0)[0] == 0)
    {
    ++numModelFaces;
    }
  //vtkIdType nextModelFace = numModelFaces + 1;
  //std::set<ModelFacetSplitRec> modelFaceSplitter;
  ModelFacetSplitMap modelFaceSplitter(numModelFaces);
  vtkIdType numCellRegions = regions.ModelRegions->GetNumberOfTuples();
  for (vtkIdType i = 0; i < numCellRegions; ++i)
    {
    vtkIdType regionIds[3];
    regions.ModelRegions->GetTupleValue(i, regionIds);
    vtkIdType currentModelFace = cell2facet(regionIds[0]);
    vtkIdType revisedModelFace = modelFaceSplitter.SplitFaceByRegion(
      currentModelFace,
      regions.Sets.Find(regionIds[1]),
      regions.Sets.Find(regionIds[2]));
    if (revisedModelFace != currentModelFace)
      {
      regionIds[0] = revisedModelFace;
      regions.ModelRegions->SetTupleValue(i, regionIds);
      modelFaceArray->SetValue(i, revisedModelFace);
      }
    }
  /*
  cout
    << numModelFaces << " faces, now " << modelFaceSplitter.NextModelFacet
    << ", " << pdIn->GetNumberOfCells() << " triangles\n";
    */
  vtkNew<vtkIdTypeArray> reconciledFaceRegionMap;
  reconciledFaceRegionMap->SetName("ModelFaceRegionsMap");
  reconciledFaceRegionMap->SetNumberOfComponents(3);
  reconciledFaceRegionMap->SetComponentName(0, "ModelFace");
  reconciledFaceRegionMap->SetComponentName(1, "BackfaceRegion");
  reconciledFaceRegionMap->SetComponentName(2, "FrontfaceRegion");
  reconciledFaceRegionMap->SetNumberOfTuples(modelFaceSplitter.NextModelFacet);
  ModelFacetSplitMap::iterator mit = modelFaceSplitter.begin();
  vtkIdType i = 0;
  for (mit = modelFaceSplitter.begin(); mit != modelFaceSplitter.end(); ++i, ++mit)
    {
    vtkIdType regionIds[3] = { mit->second, mit->first.Regions[0], mit->first.Regions[1] };
    reconciledFaceRegionMap->SetTupleValue(i, regionIds);
    }
  regions.ReconciledModelRegions = reconciledFaceRegionMap.GetPointer();
}

template<typename T, typename N>
void AssignHoles(
  vtkPolyData* pdIn, vtkPoints* holePoints,
  RegionTracker<N>& regions, T& cell2facet,
  vtkIntArray* isAHole)
{
  //cout << "Assigning holes\n";
  // For each hole point H, choose a random direction D and fire
  // a ray R=(H,D), recording intersections with each shell.
  // Find the first properly-oriented (inward-pointing) shell with
  // an odd number of hits and assign that shell to be a hole.
  //
  // NB: For safety from numerical artifacts, it would be best
  // to repeat the intersection several times for each hole point
  // and vote on the best containing shell.
  vtkNew<vtkMinimalStandardRandomSequence> rnd;
  vtkIdType numHolePts = holePoints->GetNumberOfPoints();
  for (vtkIdType i = 0; i < numHolePts; ++i)
    {
    // Get the hole point:
    vec3d h;
    holePoints->GetPoint(i, h.GetData());
    vtkIdType region = FindRegionContainingPoint(
      pdIn, h, regions, cell2facet, rnd.GetPointer());
    vtkIdType collapsedId = (region >= 0 ? regions.Collapse[region] : -1);
    if (collapsedId >= 0)
      {
      isAHole->SetValue(collapsedId, 1);
      }
    }
}

template<typename T, typename N>
void AssignRegionIDsHolesAndAttributes(
  vtkPolyData* pdIn, vtkPoints* holePtsIn, vtkPolyData* regionPtsIn,
  const std::string& vtkNotUsed(regionIdAttributeName),
  RegionTracker<N>& regions, T& cell2facet)
{
  // I. Collapse remaining union-find sets to a sequential integer numbering.
  //
  // Any pre-existing map entries will not be overwritten; we use this
  // to label the exterior shell as -1.
  //
  // Do not call regions.Sets.Find() on the new values or
  // bad things will happen quickly.
  regions.Collapse[regions.Sets.Find(regions.ExteriorId)] = -1;
  regions.Sets.CollapseIds(regions.Collapse, 0);
  regions.CollapseContainmentRelationships();

  // II. Create arrays mapping collapsed region IDs to holes and
  //     region groups.

  // Create an array to hold indexes into the user-specified groups:
  vtkNew<vtkIdTypeArray> regionGroupMembership;
  regionGroupMembership->SetName("RegionGroupMembership");
  // Don't include the exterior region:
  regionGroupMembership->SetNumberOfTuples(regions.Collapse.size() - 1);
  // Initialize regions to an invalid index (-1):
  regionGroupMembership->FillComponent(0, -1);

  // Create a matching array that specifies "hole membership"
  vtkNew<vtkIntArray> holeRegions;
  holeRegions->SetName("IsRegionAHole");
  holeRegions->SetNumberOfComponents(1);
  holeRegions->SetNumberOfTuples(regions.Collapse.size() - 1);
  holeRegions->FillComponent(0, 0);

  // III. Find which regions are marked as holes.
  if (holePtsIn && holePtsIn->GetNumberOfPoints())
    {
    AssignHoles(
      pdIn, holePtsIn, regions, cell2facet,
      holeRegions.GetPointer());
    }

  // IV. Now find which regions are in which region groups specified
  //     by the input data.
  vtkNew<vtkMinimalStandardRandomSequence> rnd;
  // Copy the region-group attribute data to our output field data:
  if (regionPtsIn)
    {
    vtkPoints* regionPts = regionPtsIn->GetPoints();
    vtkIdType numRegionPts = regionPts ? regionPts->GetNumberOfPoints() : 0;
    regions.RegionAttributes->DeepCopy(regionPtsIn->GetPointData());

    for (vtkIdType i = 0; i < numRegionPts; ++i)
      {
      // Get the region point:
      vec3d regionPt;
      regionPts->GetPoint(i, regionPt.GetData());
      // Find which region contains the point:
      vtkIdType region = FindRegionContainingPoint(
        pdIn, regionPt, regions, cell2facet, rnd.GetPointer());
      // Find the collapsed ID of the region
      vtkIdType collapsedId = (region >= 0 ? regions.Collapse[region] : -1);
      if (collapsedId >= 0)
        {
        regionGroupMembership->SetValue(collapsedId, i);
        }
      }
    }
  regions.RegionAttributes->AddArray(regionGroupMembership.GetPointer());
  regions.RegionAttributes->AddArray(holeRegions.GetPointer());

  // V. Overwrite ModelRegions with collapsed IDs.
  vtkIdType numModelFacets = regions.ModelRegions->GetNumberOfTuples();
  //cout << "TriangleRegions: (" << numModelFacets << ")\n";
  for (vtkIdType i = 0; i < numModelFacets; ++i)
    {
    vtkIdType facetRegions[3];
    regions.ModelRegions->GetTupleValue(i, facetRegions);
    for (int j = 1; j < 3; ++j)
      {
      facetRegions[j] = regions.Collapse[regions.Sets.Find(facetRegions[j])];
      }
    // Flip model-faces that have the "outside" on their backface:
    // we'll reverse the cell connectivity later in PrepareOutput.
    if (facetRegions[1] < facetRegions[2])
      {
      vtkIdType tmp = facetRegions[1];
      facetRegions[1] = facetRegions[2];
      facetRegions[2] = tmp;
      }
    regions.ModelRegions->SetTupleValue(i, facetRegions);
    //cout << i << ": " << facetRegions[0] << " " << facetRegions[1] << " " << facetRegions[2] << "\n";
    }

  // VI. Overwrite output region points with collapsed IDs.
  if (regions.RegionPoints)
    {
    vtkIdTypeArray* regionIds = vtkIdTypeArray::SafeDownCast(
      regions.RegionPoints->GetPointData()->GetArray("Region"));
    vtkIdType* reg = regionIds->GetPointer(0);
    for (vtkIdType i = 0; i < regionIds->GetNumberOfTuples(); ++i)
      {
      reg[i] = reg[i] < 0 ? reg[i] : regions.Collapse[reg[i]];
      }
    }
}

template<typename T, typename N>
void PrepareOutput(
  vtkPolyData* pdIn, RegionTracker<N>& regions, T& cell2facet,
  const char* faceGroupArrayName, vtkPolyData* pdOut, vtkTable* tabOut)
{
  pdOut->ShallowCopy(pdIn);
  vtkNew<vtkIdTypeArray> cellRegions;
  cellRegions->SetName(VTK_CELL_REGION_IDS);
  cellRegions->SetNumberOfComponents(2);
  cellRegions->SetNumberOfTuples(pdIn->GetNumberOfCells());
  vtkIdType cellIdOffset = pdIn->GetNumberOfVerts();
  cellRegions->FillComponent(0, -1);
  cellRegions->FillComponent(1, -1);
  vtkIdType polyRegions[3];
  vtkCellArray* cells =
    (N::Dimension == 2 ?
     pdIn->GetLines() : pdIn->GetPolys());
  cells->InitTraversal();
  vtkIdType npts;
  vtkIdType* conn;
  vtkIdType connOffset = 0;
  // Prepare a per-cell region array and flip cells
  // that have the "exterior" region on their backface.
  for (vtkIdType i = 0; cells->GetNextCell(npts, conn); ++i)
    {
    vtkIdType curFaceId = i + cellIdOffset;
    // Get the shell information for both co-facets of this face:
    vtkIdType modelMapEntry = regions.ModelMap[cell2facet(curFaceId)];
    regions.ModelRegions->GetTupleValue(modelMapEntry, polyRegions);
    //polyRegions[1] = regions.Sets.Find(polyRegions[1]);
    //polyRegions[2] = regions.Sets.Find(polyRegions[2]);
    //regions.ModelRegions->SetTupleValue(modelMapEntry, polyRegions);

    // Reverse the cell if its normal points inward and it is an outer shell
    // Note that this invalidates regions.EdgeNeighborhoods[...].Sense.
    if (polyRegions[1] < polyRegions[2])
      {
      cells->ReverseCell(connOffset);
      vtkIdType tmp = polyRegions[1];
      polyRegions[1] = polyRegions[2];
      polyRegions[2] = tmp;
      }

    cellRegions->SetTupleValue(curFaceId, &polyRegions[1]);
    connOffset = cells->GetTraversalLocation();
    }
  pdOut->GetCellData()->AddArray(cellRegions.GetPointer());

  // Copy the region attribute stuff to the field data.
  pdOut->GetFieldData()->ShallowCopy(regions.RegionAttributes.GetPointer());

  // Copy the model-facet boundary markers from the input if they exist:
  vtkAbstractArray* facetBoundaryMarkers =
    pdIn->GetFieldData()->GetAbstractArray(faceGroupArrayName);
  if (facetBoundaryMarkers)
    {
    pdOut->GetFieldData()->AddArray(facetBoundaryMarkers);
    }

  // If we reconciled model facets, collapse and substitute the new array
  if (
    regions.ReconciledModelRegions &&
    regions.ReconciledModelRegions->GetNumberOfTuples() > 0)
    {
    regions.ModelRegions = regions.ReconciledModelRegions;

    vtkIdType numModelFacets = regions.ModelRegions->GetNumberOfTuples();
    //cout << "ModelFaceRegions: (" << numModelFacets << ")\n";
    for (vtkIdType i = 0; i < numModelFacets; ++i)
      {
      vtkIdType facetRegions[3];
      regions.ModelRegions->GetTupleValue(i, facetRegions);
      for (int j = 1; j < 3; ++j)
        {
        facetRegions[j] = regions.Collapse[regions.Sets.Find(facetRegions[j])];
        }
      // We reversed cells above, now reverse region entries so outside is always along positive normal
      if (facetRegions[1] == -1 && facetRegions[2] != -1)
        {
        vtkIdType tmp = facetRegions[1];
        facetRegions[1] = facetRegions[2];
        facetRegions[2] = tmp;
        }
      regions.ModelRegions->SetTupleValue(i, facetRegions);
      //cout << i << ": " << facetRegions[0] << " " << facetRegions[1] << " " << facetRegions[2] << "\n";
      }
    }
  // Add the model facet -> region ID map to the field data.
  pdOut->GetFieldData()->AddArray(regions.ModelRegions.GetPointer());

  // Annotate region components
  regions.ModelRegions->SetComponentName(0, "ModelFacet");
  regions.ModelRegions->SetComponentName(1, "BackfaceRegion");
  regions.ModelRegions->SetComponentName(2, "FrontfaceRegion");
  for (int i = 0; i < 2; ++i)
    {
    cellRegions->SetComponentName(i,
      regions.ModelRegions->GetComponentName(i + 1));
    }

  tabOut->ShallowCopy(regions.ContainmentRelationships.GetPointer());
}

template<typename T, typename N>
int DiscoverRegions(
  vtkPolyData* pdIn, vtkPolyData* holesIn, vtkPolyData* regionPtsIn,
  const char* regionGroupArrayName, const char* faceGroupArrayName,
  const char* modelFaceArrayName, vtkPolyData* pdOut, vtkTable* tabOut,
  T& cell2facet, RegionTracker<N>& regions)
{
  vtkIdTypeArray* mapArray = NULL;

  InitializeRegions(pdIn, regions, cell2facet);
  if (pdIn->GetNumberOfCells())
    {
    // Merge regions with adjacent (dim-2)-boundaries.
    MergeRegions(pdIn, regions, cell2facet);

    // Find which regions (if any) are contained in other regions
    // Merge them as appropriate.
    DiscoverNestings(pdIn, regions, cell2facet);

    // It is possible that some model faces are "oversubscribed"
    // (i.e., they are bordant to more than 2 regions). This can
    // only occur if model faces were not initially simplices
    // but were provided to us after tessellation -- defined by
    // way of a cell-data scalar specifying the original model
    // face.
    // This is corrected by splitting (d-1)-simplices in the
    // tessellation of model faces into groups that share the
    // same regions.
    //
    // This will alter the array specified by modelFaceArrayName
    // by changing existing model face assignments as required.
    // It will also change entries in regions.ModelRegions.
    mapArray =
      modelFaceArrayName ?
      vtkIdTypeArray::SafeDownCast(
        pdIn->GetCellData()->GetArray(modelFaceArrayName)) :
      NULL;
    if (mapArray)
      {
      vtkDiscoverRegionsFacetFromCellArray<vtkIdTypeArray> reconciler;
      reconciler.Array = mapArray;
      ReconcileModelFacets(pdIn, modelFaceArrayName, regions, reconciler);
      }
    }

  // See if we have enough information to mark regions
  // as belonging to domain sets.
  vtkPoints* regionPts = regionPtsIn ? regionPtsIn->GetPoints() : NULL;
  if (
    !regionGroupArrayName || !regionGroupArrayName[0] ||
    !regionPts || !regionPts->GetNumberOfPoints())
    {
    regionPtsIn = NULL;
    }

  // Points identifying regions as holes.
  vtkPoints* holePointsIn = holesIn ? holesIn->GetPoints() : NULL;

  // Map region IDs into a 0-based integer sequence.
  AssignRegionIDsHolesAndAttributes(
    pdIn, holePointsIn, regionPtsIn, regionGroupArrayName ? regionGroupArrayName : "", regions, cell2facet);

  // Copy region IDs to output dataset
  PrepareOutput(pdIn, regions, cell2facet, faceGroupArrayName, pdOut, tabOut);
  return 1;
}

int vtkDiscoverRegions::FillInputPortInformation(
  int port,
  vtkInformation* info)
{
  if (!this->Superclass::FillInputPortInformation(port, info))
    {
    return 0;
    }
  if (port > 0)
    {
    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
    }
  return 1;
}

int vtkDiscoverRegions::FillOutputPortInformation(
  int port,
  vtkInformation* info)
{
  if (port == 1)
    {
    info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable");
    }
  else if (!this->Superclass::FillOutputPortInformation(port, info))
    {
    return 0;
    }
  return 1;
}

int vtkDiscoverRegions::RequestData(
    vtkInformation* vtkNotUsed(req),
    vtkInformationVector** inInfo,
    vtkInformationVector* outInfo)
{
  // get the input and output
  vtkPolyData* pdIn = vtkPolyData::GetData(inInfo[0], 0);
  vtkPolyData* holesIn = vtkPolyData::GetData(inInfo[1], 0);
  vtkPolyData* regionPtsIn = vtkPolyData::GetData(inInfo[2], 0);
  vtkPolyData* pdOut = vtkPolyData::GetData(outInfo, 0);
  vtkTable* tabOut = vtkTable::GetData(outInfo, 1);
  vtkPolyData* regionPtsOut =
    this->GenerateRegionInteriorPoints ?
    vtkPolyData::GetData(outInfo, 2) : NULL;

  /*
  vtkPolyData* pdIn = this->GetPolyDataInput(0);
  vtkPolyData* holesIn = this->GetPolyDataInput(1);
  vtkPolyData* regionPtsIn = this->GetPolyDataInput(2);
  vtkPolyData* pdOut = this->GetOutput();
  */

  vtkIdTypeArray* mapArray =
    this->ModelFaceArrayName && !this->ReportRegionsByModelFace ?
      vtkIdTypeArray::SafeDownCast(
        pdIn->GetCellData()->GetArray(this->ModelFaceArrayName)) :
      NULL;

  int dimension = -1;
  if (!pdIn->GetNumberOfPolys() && !pdIn->GetNumberOfStrips())
    {
    if (pdIn->GetNumberOfLines() > 0)
      {
      dimension = 2;
      }
    }
  else
    {
    dimension = 3;
    }

  int status;
  switch (dimension)
    {
  case 2:
      {
      RegionTracker<Vertexhoods> regionTracker;
      regionTracker.RegionPoints = regionPtsOut;
      if (mapArray)
        {
        vtkDiscoverRegionsFacetFromCellArray<vtkIdTypeArray> cell2facet;
        cell2facet.Array = mapArray;
        status = DiscoverRegions(
          pdIn, holesIn, regionPtsIn,
          this->RegionGroupArrayName,
          this->FaceGroupArrayName,
          this->ModelFaceArrayName,
          pdOut, tabOut, cell2facet, regionTracker);
        }
      else
        { // NULL or invalid name: one facet per poly
        vtkDiscoverRegionsOneFacetPerCell cell2facet;
        status = DiscoverRegions(
          pdIn, holesIn, regionPtsIn,
          this->RegionGroupArrayName,
          this->FaceGroupArrayName,
          this->ModelFaceArrayName,
          pdOut, tabOut, cell2facet, regionTracker);
        }
      }
    break;
  case 3:
      {
      RegionTracker<Edgerhoods> regionTracker;
      regionTracker.RegionPoints = regionPtsOut;
      if (mapArray)
        {
        vtkDiscoverRegionsFacetFromCellArray<vtkIdTypeArray> cell2facet;
        cell2facet.Array = mapArray;
        status = DiscoverRegions(
          pdIn, holesIn, regionPtsIn,
          this->RegionGroupArrayName,
          this->FaceGroupArrayName,
          this->ModelFaceArrayName,
          pdOut, tabOut, cell2facet, regionTracker);
        }
      else
        { // NULL or invalid name: one facet per poly
        vtkDiscoverRegionsOneFacetPerCell cell2facet;
        status = DiscoverRegions(
          pdIn, holesIn, regionPtsIn,
          this->RegionGroupArrayName,
          this->FaceGroupArrayName,
          this->ModelFaceArrayName,
          pdOut, tabOut, cell2facet, regionTracker);
        }
      }
    break;
  default:
    //vtkErrorMacro("Unable to determine problem dimension.");
    return 1;
    }

  return 1;
}
