vtkVoronoiTile.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 vtkVoronoiTile_h
#define vtkVoronoiTile_h
 
#include "vtkDoubleArray.h"          // Support double points
#include "vtkFiltersMeshingModule.h" // For export macro
#include "vtkMath.h"                 // Euclidean norm
#include "vtkVoronoiCore.h"          // Enum: clip intersection status
 
#include <vector> // array support
 
VTK_ABI_NAMESPACE_BEGIN
 
class vtkPolyData;
class vtkSpheres;
 
//======= Define the convex polygon class used to produce Voronoi tiles.
 
// The data structure for representing a Voronoi tile vertex and implicitly,
// the connected Voronoi tile edge. The tile vertex has a position X, and the
// current value of the half-space clipping function. In the counterclockwise
// direction, the NeiId refers to the point id in the neighboring tile that,
// together with this tile's point id, produced a tile edge between the two
// points (i.e., a spoke).
struct vtkTilePoint
{
  double X[2];     // position of this vertex
  vtkIdType NeiId; // generating point id for the associated edge
  double Val;      // current value of the current half-space clipping function
  double R2;       // Radius**2 of circumcircle / flower petal
 
  vtkTilePoint(double tileX[2], double x[2], vtkIdType neiId)
    : X{ x[0], x[1] }
    , NeiId(neiId)
    , Val(0.0)
  {
    this->R2 = vtkMath::Distance2BetweenPoints2D(X, tileX);
  }
 
  vtkTilePoint(const vtkTilePoint& v) = default;
};
 
// Typedefs defined for convenience.
using PointRingType = std::vector<vtkTilePoint>;
using PointRingIterator = std::vector<vtkTilePoint>::iterator;

class VTKFILTERSMESHING_EXPORT vtkVoronoiTile
{
public:
  vtkVoronoiTile()
    : PtId(-1)
    , X{ 0, 0, 0 } // z-component specifies location in z-plane
    , NumClips(0)
    , PruneTolerance(1.0e-13)
    , RecomputeCircumFlower(true)
    , RecomputePetals(true)
    , CircumFlower2(0.0)
    , MinRadius2(0.0)
    , MaxRadius2(0.0)
  {
    // Preallocate some space
    this->Points.reserve(256);
    this->NewPoints.reserve(256);
 
    // Supporting data structures
    this->SortP.reserve(256);
    this->Petals = vtkSmartPointer<vtkDoubleArray>::New();
    this->Petals->SetNumberOfComponents(3); // x-y-R2
    this->Petals->Allocate(256);            // initial allocation
  }

  void Initialize(vtkIdType genPtId, const double genPt[3], double bds[4]);

  void Initialize(
    vtkIdType ptId, const double x[3], vtkPoints* pts, vtkIdType nPts, const vtkIdType* p);

  ClipIntersectionStatus Clip(vtkIdType neiPtId, const double neiPt[2]);

  double GetCircumFlower2()
  {
    if (this->RecomputeCircumFlower)
    {
      this->ComputeCircumFlower();
    }
    return this->CircumFlower2;
  }
  bool InCircumFlower(double r2) // radius**2 of point from generator
  {
    // Only recompute the circumflower if necessary; that is, when
    // a maximal point is eliminated by a plane clip.
    if (this->RecomputeCircumFlower)
    {
      this->ComputeCircumFlower();
    }
    return (r2 <= this->CircumFlower2);
  }
  bool InFlower(const double x[2]);
  void UpdatePetals(double cf2);
  vtkDoubleArray* GetPetals()
  {
    if (this->RecomputePetals)
    {
      this->UpdatePetals(this->CircumFlower2);
    }
    return (this->Petals->GetNumberOfTuples() > 0 ? this->Petals : nullptr);
  }

  void ProducePolyData(vtkPolyData* pd, vtkSpheres* spheres);
  void ProducePolyData(vtkPolyData* pd) { this->ProducePolyData(pd, nullptr); }

  vtkIdType GetGeneratorPointId() { return this->PtId; }
  double* GetGeneratorPosition() { return this->X; }
  vtkIdType GetNumberOfPoints() { return this->Points.size(); }
  const PointRingType& GetPoints() { return this->Points; }
 
  // Data members purposely left public for using classes to extract
  // information.
 
  // Information used to define the polyhedron- its generating point id and
  // position, plus region classification. Indicate whether degenerate faces
  // (i.e., those having ~zero area) can be deleted (i.e., pruned).
  vtkIdType PtId;        // Generating point id
  double X[3];           // Generating point position. Note X[2] is z-plane position
  vtkIdType NumClips;    // The total number of clip operations since Initialize()
  double PruneTolerance; // Specify the prune tolerance
 
  // These data members represent the constructed polygon.
  PointRingType Points;    // counterclockwise ordered loop of points/vertices defining the tile
  PointRingType NewPoints; // accumulate new points/vertices to construct the tile
 
protected:
  // Convenience method for moving around the modulo ring of the tile
  // vertices.
  PointRingIterator Next(PointRingIterator& itr)
  {
    if (itr == (this->Points.end() - 1))
    {
      return this->Points.begin();
    }
    return (itr + 1);
  }
 
  // The core geometric intersection operation.
  ClipIntersectionStatus IntersectWithLine(double origin[2], double normal[2], vtkIdType neiPtId);
 
  // These tolerances are for managing degeneracies.
  double Tol;
  double Tol2;
 
  // Indicate whether the Voronoi circumflower needs recomputing, and
  // keep track of the current circumflower and related information.
  void ComputeCircumFlower();
  bool RecomputeCircumFlower;
  bool RecomputePetals;
  double CircumFlower2;
  double MinRadius2;
  double MaxRadius2;
  std::vector<vtkTilePoint*> SortP;       // Points sorted on radius**2
  vtkSmartPointer<vtkDoubleArray> Petals; // Flower petals w/ radii > annular radius
 
}; // vtkVoronoiTile
 
//------------------------------------------------------------------------------
inline void vtkVoronoiTile::ComputeCircumFlower()
{
  // Compute the circumflower, and compute some info about
  // the flower radii.
  this->MinRadius2 = VTK_FLOAT_MAX;
  this->MaxRadius2 = VTK_FLOAT_MIN;
 
  // Determine the circumflower and minimal sphere radius by
  // checking against each of the flower petals.
  for (const auto& p : this->Points)
  {
    double r2 = p.R2;
    this->MinRadius2 = std::min(r2, this->MinRadius2);
    this->MaxRadius2 = std::max(r2, this->MaxRadius2);
  }
  this->CircumFlower2 = (4.0 * this->MaxRadius2); // (2*(max petal radius))**2
  this->RecomputeCircumFlower = false;            // circumflower is up to date
}
 
//------------------------------------------------------------------------------
inline bool vtkVoronoiTile::InFlower(const double x[2])
{
  // Check against the flower petals
  for (const auto& p : this->Points)
  {
    double r2 = vtkMath::Distance2BetweenPoints2D(p.X, x);
    if (r2 <= p.R2)
    {
      return true;
    }
  }
 
  // Point not in the flower
  return false;
}
 
VTK_ABI_NAMESPACE_END
#endif
// VTK-HeaderTest-Exclude: vtkVoronoiTile.h