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->ReserveValues(256);       // initial allocation
  }

  vtkVoronoiTile(const 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);

    // Each tile owns its own vtkDoubleArray.
    this->SortP.reserve(256);
    this->Petals = vtkSmartPointer<vtkDoubleArray>::New();
    this->Petals->SetNumberOfComponents(3); // x-y-R2
    this->Petals->ReserveTuples(256);       // initial allocation
  }

  vtkVoronoiTile& operator=(const vtkVoronoiTile&) { return *this; }

  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