Method ComputeBoundingSphere

Common/DataModel/vtkCell.cxx:157–285

------------------------------------------------------------------------------

//------------------------------------------------------------------------------
double vtkCell::ComputeBoundingSphere(double center[3]) const
{
  // We do easy cases first for number of points <= 4
  switch (this->Points->GetNumberOfPoints())
  {
    case 0:
      center[0] = std::numeric_limits<double>::quiet_NaN();
      center[1] = std::numeric_limits<double>::quiet_NaN();
      center[2] = std::numeric_limits<double>::quiet_NaN();
      return std::numeric_limits<double>::quiet_NaN();
    case 1:
      this->Points->GetPoint(0, center);
      return 0.0;
    case 2:
    {
      auto points = vtk::DataArrayTupleRange(this->Points->GetData());
      auto p0 = points[0], p1 = points[1];
      center[0] = 0.5 * (p0[0] + p1[0]);
      center[1] = 0.5 * (p0[1] + p1[1]);
      center[2] = 0.5 * (p0[2] + p1[2]);
      return vtkMath::Distance2BetweenPoints(center, p0);
    }
    case 3:
    {
      if (!vtkTriangle::ComputeCentroid(this->Points, nullptr, center))
      {
        break;
      }
      auto points = vtk::DataArrayTupleRange(this->Points->GetData());
      return vtkMath::Distance2BetweenPoints(center, points[0]);
    }
    case 4:
    {
      if (!vtkTetra::ComputeCentroid(this->Points, nullptr, center))
      {
        break;
      }
      auto points = vtk::DataArrayTupleRange(this->Points->GetData());
      return vtkMath::Distance2BetweenPoints(center, points[0]);
    }
    default:
      break;
  }

  // For more complex cells, we follow Ritter's bounding sphere algorithm
  // 1. Pick a point x (first point in our case) in the cell, and look for
  //    a point y the furthest from x.
  // 2. Look for a point z the furthest from y.
  // 3. Create a sphere centered at [z,y] with appropriate radius
  // 4. Until all points are not in the sphere, take a point outside the sphere,
  //    and update the sphere to include former sphere + this point

  auto points = vtk::DataArrayTupleRange(this->Points->GetData());
  using ConstRefType = typename decltype(points)::ConstTupleReferenceType;

  ConstRefType x = points[0];
  vtkIdType yid = 1, zid = 0;