Method IntersectWithLine

Common/DataModel/vtkCellLocator.cxx:135–316

------------------------------------------------------------------------------ Return intersection point (if any) AND the cell which was intersected by finite line.

// Return intersection point (if any) AND the cell which was intersected by
// finite line.
int vtkCellLocator::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
  double x[3], double pcoords[3], int& subId, vtkIdType& cellId, vtkGenericCell* cell)
{
  this->BuildLocator();
  if (this->Tree == nullptr)
  {
    return 0;
  }
  double* bounds = this->Bounds;
  double* h = this->H;
  double t0, t1, x0[3], x1[3], tHitCell;
  double hitCellBoundsPosition[3], cellBounds[6], *cellBoundsPtr;
  double octantBounds[6];
  int prod = this->NumberOfDivisions * this->NumberOfDivisions;
  vtkIdType cellIdBest = -1, cId, i, idx, numberOfCellsInBucket;
  int ijk[3], ijkEnd[3];
  int plane0, plane1, subIdBest = -1, hitCellBounds;
  double tBest = VTK_FLOAT_MAX, xBest[3], pCoordsBest[3], step[3], next[3], tMax[3], tDelta[3];
  double rayDir[3];
  vtkMath::Subtract(p2, p1, rayDir);
  int leafStart = this->NumberOfOctants -
    this->NumberOfDivisions * this->NumberOfDivisions * this->NumberOfDivisions;

  // Make sure the bounding box of the locator is hit. Also, determine the
  // entry and exit points into and out of the locator. This is used to
  // determine the bins where the ray starts and ends.
  cellId = -1;
  subId = 0;
  if (vtkBox::IntersectWithLine(bounds, p1, p2, t0, t1, x0, x1, plane0, plane1) == 0)
  {
    return 0; // No intersections possible, line is outside the locator
  }

  // Initialize intersection query array if necessary. This is done
  // locally to ensure thread safety.
  std::vector<bool> cellHasBeenVisited(this->DataSet->GetNumberOfCells(), false);

  // Get the i-j-k point of intersection and bin index. This is
  // clamped to the boundary of the locator.
  this->GetBucketIndices(x0, ijk);
  this->GetBucketIndices(x1, ijkEnd);
  idx = leafStart + ijk[0] + ijk[1] * this->NumberOfDivisions + ijk[2] * prod;

  // Set up some traversal parameters for traversing through bins
  step[0] = (rayDir[0] >= 0.0) ? 1.0 : -1.0;
  step[1] = (rayDir[1] >= 0.0) ? 1.0 : -1.0;
  step[2] = (rayDir[2] >= 0.0) ? 1.0 : -1.0;

  // If the ray is going in the negative direction, then the next voxel boundary
  // is on the "-" direction so we stay in the current voxel.
  next[0] = bounds[0] + h[0] * (rayDir[0] >= 0.0 ? (ijk[0] + step[0]) : ijk[0]);
  next[1] = bounds[2] + h[1] * (rayDir[1] >= 0.0 ? (ijk[1] + step[1]) : ijk[1]);
  next[2] = bounds[4] + h[2] * (rayDir[2] >= 0.0 ? (ijk[2] + step[2]) : ijk[2]);

  tMax[0] = (rayDir[0] != 0.0) ? (next[0] - x0[0]) / rayDir[0] : VTK_FLOAT_MAX;
  tMax[1] = (rayDir[1] != 0.0) ? (next[1] - x0[1]) / rayDir[1] : VTK_FLOAT_MAX;
  tMax[2] = (rayDir[2] != 0.0) ? (next[2] - x0[2]) / rayDir[2] : VTK_FLOAT_MAX;