vtkVoronoiCore.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#ifndef vtkVoronoiCore_h
#define vtkVoronoiCore_h

#include "vtkAlgorithm.h" // check abort status if embedded in filter
#include "vtkIntArray.h"  // for representing segmented region ids
#include "vtkSMPTools.h"  // SMP parallel processing

#include <algorithm> // for std::sort
#include <iostream>
#include <random> // random generation of colors
#include <vector>

VTK_ABI_NAMESPACE_BEGIN

enum class ClipIntersectionStatus
{
  NoIntersection = 0,
  Intersection = 1,
  Pruned = 2,
  Numeric = 3
};

enum vtkSpokeClassification
{
  BACKWARD_SPOKE = 0,  // Bit 0: Backward spoke
  FORWARD_SPOKE = 1,   // Bit 0: Forward spoke
  REGION_BOUNDARY = 2, // Bit 1: Region boundary spoke
  DOMAIN_BOUNDARY = 4, // Bit 2: Domain boundary spoke
  PRUNED = 8,          // Bit 3: Spoke is pruned (deleted)
};

struct vtkVoronoiSpoke
{
  vtkIdType NeiId;              // Id of the wheel that the spoke is connected to (wheelId,NeiId)
  unsigned char Classification; // Indicate the classification of this spoke
  vtkVoronoiSpoke()
    : NeiId(-1)
    , Classification(0)
  {
  }
  vtkVoronoiSpoke(vtkIdType neiId, unsigned char classification)
    : NeiId(neiId)
    , Classification(classification)
  {
  }
}; // vtkVoronoiSpoke

using vtkVoronoiSpokesType = std::vector<vtkVoronoiSpoke>;
using vtkVoronoiSpokesIterator = vtkVoronoiSpokesType::iterator;
using vtkVoronoiWheelsType = std::vector<vtkIdType>;

// Gather spokes into a wheel. Define some basic operators.  Note that every
// wheel is associated with an input (tile/hull generating) point. So access
// to the wheel and its associated spokes is via point id.
struct vtkVoronoiWheel
{
  vtkVoronoiWheelsType& Wheels; // The composited array of wheels
  vtkVoronoiSpokesType& Spokes; // The composited array of spokes
  vtkIdType Id;                 // The associated point/tile id: with wheelId == pointId
  int NumSpokes;                // The number of emanating spokes
  vtkVoronoiSpoke* LocalSpokes; // The array of spokes connected to this wheel

  // Default instantiation.
  vtkVoronoiWheel(vtkVoronoiWheelsType& wheels, vtkVoronoiSpokesType& spokes)
    : Wheels(wheels)
    , Spokes(spokes)
    , Id(-1)
    , NumSpokes(0)
    , LocalSpokes(nullptr)
  {
  }

  // Setup the wheel for queries: an efficient form that does not require
  // repeated wheel instantiation.
  vtkVoronoiSpoke* Initialize(vtkIdType id, int& numSpokes)
  {
    this->Id = id;
    this->NumSpokes = numSpokes = (this->Wheels[id + 1] - this->Wheels[id]);
    this->LocalSpokes =
      ((size_t)this->Wheels[id] < this->Spokes.size() ? &(this->Spokes[this->Wheels[id]])
                                                      : nullptr);
    return this->LocalSpokes;
  }
}; // vtkVoronoiWheel

struct vtkVoronoiAdjacencyGraph
{
  vtkVoronoiWheelsType Wheels; // Wheel/spokes data structure: offset array into spokes
  vtkVoronoiSpokesType Spokes; // Spokes / edges with classification

  void Initialize(vtkIdType numWheels, vtkIdType numSpokes);
  vtkVoronoiWheelsType& GetWheels() { return this->Wheels; }
  vtkVoronoiSpokesType& GetSpokes() { return this->Spokes; }
  vtkIdType GetNumberOfWheels() { return (this->Wheels.size() - 1); }
  vtkIdType GetNumberOfSpokes() { return this->Spokes.size(); }
  bool IsSpoke(vtkIdType pt0, vtkIdType pt1);
  vtkIdType GetNumberOfSpokes(vtkIdType ptId); // #spokes for a specified wheel
  vtkVoronoiSpoke* GetSpokes(vtkIdType ptId, vtkIdType& numSpokes);
  vtkIdType GetWheelOffset(vtkIdType ptId) { return this->Wheels[ptId]; }
  static void CountFaces(const vtkVoronoiSpoke* spokes, int numSpokes, int& numDomainBoundaryFaces,
    int& numRegionBoundaryFaces, int& numForwardFaces);

  bool Validate();

  struct ValidateAdjacencyGraph
  {
    vtkVoronoiAdjacencyGraph& Graph;
    vtkIdType NumInvalid;

    ValidateAdjacencyGraph(vtkVoronoiAdjacencyGraph& graph)
      : Graph(graph)
      , NumInvalid(0)
    {
    }

    // Keep track whether threads are non-degenerate.
    vtkSMPThreadLocal<vtkIdType> ThreadNumInvalid;

    // vtkSMPTools threaded interface
    void Initialize();
    void operator()(vtkIdType wheelId, vtkIdType endWheelId);
    void Reduce();
  }; // ValidateAdjacencyGraph
}; // vtkVoronoiAdjacencyGraph

struct vtkVoronoiHullVertex // 3D hull vertices
{
  double X[3];
  vtkVoronoiHullVertex(double x, double y, double z)
    : X{ x, y, z }
  {
  }
  vtkVoronoiHullVertex(const double x[3])
    : X{ x[0], x[1], x[2] }
  {
  }
};
using vtkVoronoiHullVertexType = std::vector<vtkVoronoiHullVertex>;

struct vtkVoronoiTileVertex // 2D tile vertices
{
  double X[2];
  vtkVoronoiTileVertex(double x, double y)
    : X{ x, y }
  {
  }
  vtkVoronoiTileVertex(const double x[2])
    : X{ x[0], x[1] }
  {
  }
};
using vtkVoronoiTileVertexType = std::vector<vtkVoronoiTileVertex>;

struct vtkVoronoiTopoCoord3D
{
  std::array<vtkIdType, 4> Ids;

  vtkVoronoiTopoCoord3D()
    : Ids{ 0 }
  {
  }
  vtkVoronoiTopoCoord3D(vtkIdType p0, vtkIdType p1, vtkIdType p2, vtkIdType ptId)
    : Ids{ p0, p1, p2, ptId }
  {
    std::sort(this->Ids.data(), this->Ids.data() + this->Ids.size());
  }
  vtkVoronoiTopoCoord3D(const vtkVoronoiTopoCoord3D& tt) { this->Ids = tt.Ids; }

  bool operator<(const vtkVoronoiTopoCoord3D& tuple) const { return this->Ids < tuple.Ids; }
};
using vtkVoronoiTopoCoords3DType = std::vector<vtkVoronoiTopoCoord3D>;

struct vtkVoronoiTopoCoord2D
{
  std::array<vtkIdType, 3> Ids;

  vtkVoronoiTopoCoord2D()
    : Ids{ 0 }
  {
  }
  vtkVoronoiTopoCoord2D(vtkIdType p0, vtkIdType p1, vtkIdType ptId)
    : Ids{ p0, p1, ptId }
  {
    std::sort(this->Ids.data(), this->Ids.data() + this->Ids.size());
  }
  vtkVoronoiTopoCoord2D(const vtkVoronoiTopoCoord2D& tt) { this->Ids = tt.Ids; }

  bool operator<(const vtkVoronoiTopoCoord2D& tuple) const { return this->Ids < tuple.Ids; }
};
using vtkVoronoiTopoCoords2DType = std::vector<vtkVoronoiTopoCoord2D>;

struct vtkVoronoiMergeTuple3D : public vtkVoronoiTopoCoord3D
{
  vtkIdType PtId; // the id of the hull vertex

  vtkVoronoiMergeTuple3D()
    : PtId(-1)
  {
  }
  bool operator!=(const vtkVoronoiMergeTuple3D& mt) const { return (this->Ids != mt.Ids); }
}; // vtkVoronoiMergeTuple3D

struct vtkVoronoiMergeTuple2D : public vtkVoronoiTopoCoord2D
{
  vtkIdType PtId; // the id of the tile vertex

  vtkVoronoiMergeTuple2D()
    : PtId(-1)
  {
  }
  bool operator!=(const vtkVoronoiMergeTuple2D& mt) const { return (this->Ids != mt.Ids); }
}; // vtkVoronoiMergeTuple2D

using vtkMergeTupleOffsets = std::vector<vtkIdType>; // offsets into merged tuples
using vtkMergeTuples3DType = std::vector<vtkVoronoiMergeTuple3D>;
using vtkMergeTuples2DType = std::vector<vtkVoronoiMergeTuple2D>;

using vtkMergeMapType = std::vector<vtkIdType>;

using vtkVoronoiCellConnType = std::vector<vtkIdType>;

struct vtkVoronoiBatchManager
{
  vtkIdType Num;        // Number of total items (e.g., points) to process
  vtkIdType BatchSize;  // The desired batch size (clamped by Num)
  vtkIdType NumBatches; // The total number of batches to process
  vtkVoronoiBatchManager(vtkIdType num, vtkIdType batchSize)
    : Num(num)
    , BatchSize(batchSize)
  {
    this->NumBatches = static_cast<vtkIdType>(ceil(static_cast<double>(num) / batchSize));
  }
  vtkIdType GetNumberOfBatches() const { return this->NumBatches; }
  vtkIdType GetBatchItemRange(vtkIdType batchNum, vtkIdType& startId, vtkIdType& endId) const
  {
    startId = batchNum * this->BatchSize;
    endId = startId + this->BatchSize;
    endId = (endId > this->Num ? this->Num : endId);
    return (endId - startId);
  }
}; // vtkVoronoiBatchManager

using vtkBatchIdsType = std::vector<vtkIdType>;

// Convenience function: convert input labels/region ids/scalars to signed int.
// The Voronoi classes expect signed int region labels.
inline vtkSmartPointer<vtkIntArray> ConvertRegionLabels(vtkDataArray* inScalars)
{
  vtkNew<vtkIntArray> rIds;
  rIds->SetNumberOfTuples(inScalars->GetNumberOfTuples());
  rIds->DeepCopy(inScalars);

  return rIds;
}

struct vtkVoronoiAbortCheck
{
  vtkAlgorithm* Filter;
  bool IsFirst;
  vtkIdType CheckAbortInterval;

  vtkVoronoiAbortCheck(vtkIdType start, vtkIdType end, vtkAlgorithm* filter)
    : Filter(filter)
  {
    this->IsFirst = vtkSMPTools::GetSingleThread();
    this->CheckAbortInterval = std::min((end - start) / 10 + 1, (vtkIdType)1000);
  }

  bool operator()(vtkIdType id)
  {
    if (this->Filter && this->IsFirst && !(id % this->CheckAbortInterval))
    {
      this->Filter->CheckAbort();
      return (this->Filter->GetAbortOutput() ? true : false);
    }
    return false;
  }
}; // vtkVoronoiAbortCheck

// Use system <random> - create a simple convenience class. This generates
// random color indices [0,64).
struct vtkVoronoiRandomColors
{
  std::mt19937 RNG;
  std::uniform_int_distribution<int> Dist;
  vtkVoronoiRandomColors() { this->Dist.param(typename decltype(this->Dist)::param_type(0, 64)); }
  void Seed(vtkIdType s) { this->RNG.seed(s); }
  vtkIdType Next() { return this->Dist(RNG); }
};

// Use system <random> - create a simple convenience class. This generates
// random real values [0,1).
struct vtkVoronoiRandom01Range
{
  std::mt19937 RNG;
  std::uniform_real_distribution<double> Dist;
  vtkVoronoiRandom01Range()
  {
    this->Dist.param(typename decltype(this->Dist)::param_type(0.0, 1.0));
  }
  void Seed(vtkIdType s) { this->RNG.seed(s); }
  double Next() { return this->Dist(RNG); }
};

// A convenience class and methods to randomly perturb (alternatively:
// joggle, jitter, or jiggle) point positions. Such jittering (even if very
// small) significantly improves the numerical stability of Voronoi and
// Delaunay computations. Terminology note: while joggle is not a widely used
// term (typically in computer graphics the words perturb, jitter, and jiggle
// are found more frequently), it is used by the popular QuickHull (and
// associated implementation QHull qhull.org)--so the terminology was
// adopted by this work. Implementation note: these methods might be better
// added to vtkMath since they can be used by other classes. Once they are
// demonstrated to be stable, they may be moved.
struct vtkVoronoiJoggle
{
  // Joggle a single point at input position xIn to produce the output position
  // xOut (xIn and xOut may be computed in-place). The radius is the allowable
  // range of joggle in the sphere. A sequence is provided, assumed properly
  // initialized, to produce random (0,1) values. Note that if this method is
  // invoked in a thread, separate sequence instantiations (one per thread)
  // should be provided.
  static void JoggleXYZ(
    double xIn[3], double xOut[3], double radius, vtkVoronoiRandom01Range& sequence)
  {
    double cosphi = 1 - 2 * sequence.Next();
    double sinphi = sqrt(1 - cosphi * cosphi);
    double rho = radius * pow(sequence.Next(), 0.33333333);
    double R = rho * sinphi;
    double theta = 2.0 * vtkMath::Pi() * sequence.Next();
    xOut[0] = xIn[0] + R * cos(theta);
    xOut[1] = xIn[1] + R * sin(theta);
    xOut[2] = xIn[2] + rho * cosphi;
  }

  // Joggle a single point at input position xIn to produce the output
  // position xOut (xIn and xOut may be computed in-place). The radius is the
  // allowable range of joggle in the circle in the x-y plane. A sequence is
  // provided, assumed properly initialized, to produce random (0,1) values.
  // Note that if this method is invoked in a thread, separate sequence
  // instantiations (one per thread) should be provided.
  static void JoggleXY(
    double xIn[3], double xOut[3], double radius, vtkVoronoiRandom01Range& sequence)
  {
    double R = radius * sequence.Next();
    double theta = 2.0 * vtkMath::Pi() * sequence.Next();
    xOut[0] = xIn[0] + R * cos(theta);
    xOut[1] = xIn[1] + R * sin(theta);
    xOut[2] = xIn[2];
  }

  // Joggle a single point at input position xIn to produce the output
  // position xOut (xIn and xOut may be computed in-place). The radius is the
  // allowable range of joggle in the circle in the x-z plane. A sequence is
  // provided, assumed properly initialized, to produce random (0,1) values.
  // Note that if this method is invoked in a thread, separate sequence
  // instantiations (one per thread) should be provided.
  static void JoggleXZ(
    double xIn[3], double xOut[3], double radius, vtkVoronoiRandom01Range& sequence)
  {
    double R = radius * sequence.Next();
    double theta = 2.0 * vtkMath::Pi() * sequence.Next();
    xOut[0] = xIn[0] + R * cos(theta);
    xOut[1] = xIn[1];
    xOut[2] = xIn[2] + R * sin(theta);
  }

  // Joggle a single point at input position xIn to produce the output
  // position xOut (xIn and xOut may be computed in-place). The radius is the
  // allowable range of joggle in the circle in the y-z plane. A sequence is
  // provided, assumed properly initialized, to produce random (0,1) values.
  // Note that if this method is invoked in a thread, separate sequence
  // instantiations (one per thread) should be provided.
  static void JoggleYZ(
    double xIn[3], double xOut[3], double radius, vtkVoronoiRandom01Range& sequence)
  {
    double R = radius * sequence.Next();
    double theta = 2.0 * vtkMath::Pi() * sequence.Next();
    xOut[0] = xIn[0];
    xOut[1] = xIn[1] + R * cos(theta);
    xOut[2] = xIn[2] + R * sin(theta);
  }

  // Joggle a 3D vector. Specify an input normalized vector (vecIn), a angle
  // of rotation theta (in radians, 0<theta<Pi/2), and a random
  // sequence. Produce on output a normalized vector (vecOut) -- vecIn and
  // vecOut may be computed in place. Note that if this method is invoked in
  // a thread, separate sequence instantiations (one per thread) should be
  // provided. The approach used here is to define a disk located at the tip
  // of the (normalized) vecIn (the disc is normal to vecIn), and randomly
  // select a vector from the edge/along the perimeter of the disc.
  static void JoggleNormal(
    double vecIn[3], double vecOut[3], double theta, vtkVoronoiRandom01Range& sequence)
  {
    // Find a vector orthogonal to the input vector.
    double perp[3];
    int iMax = (std::fabs(vecIn[0]) > std::fabs(vecIn[1]) ? 0 : 1);
    iMax = (std::fabs(vecIn[iMax]) > std::fabs(vecIn[2]) ? iMax : 2);
    perp[(iMax + 1) % 3] = 1.0;
    perp[(iMax + 2) % 3] = 0.0;
    perp[iMax] = -(vecIn[(iMax + 1) % 3] / vecIn[iMax]);
    vtkMath::Normalize(perp);

    // Now find third orthogonal vector. Since vecIn and perp are normalized and
    // orthogonal, the resultant cross also a unit vector.
    double cross[3];
    vtkMath::Cross(vecIn, perp, cross);

    // Compute a random vector
    double radius = theta / (2.0 * vtkMath::Pi()); // small angle approximation
    double t = 2.0 * vtkMath::Pi() * sequence.Next();
    double x = radius * cos(t);
    double y = radius * sin(t);

    // Construct a vector along the disc/cone edge
    vecOut[0] = vecIn[0] + x * perp[0] + y * cross[0];
    vecOut[1] = vecIn[1] + x * perp[1] + y * cross[1];
    vecOut[2] = vecIn[2] + x * perp[2] + y * cross[2];
    vtkMath::Normalize(vecOut);
  }
}; // vtkVoronoiJoggle

VTK_ABI_NAMESPACE_END
#include "vtkVoronoiCore.txx"

#endif
// VTK-HeaderTest-Exclude: vtkVoronoiCore.h