Method IntersectWithLine

Common/DataModel/vtkStaticCellLocator.cxx:598–787

// the IntersectWithLine method for more information on voxel traversal.
template <typename T>
int CellProcessor<T>::IntersectWithLine(const double p1[3], const double p2[3], double tol,
  vtkPoints* points, vtkIdList* cellIds, vtkGenericCell* cell)
{
  // Initialize the list of points/cells
  if (points)
  {
    points->Reset();
  }
  if (cellIds)
  {
    cellIds->Reset();
  }

  // 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.
  double* bounds = this->Binner->Bounds;
  double t0, t1, x0[3], x1[3], rayDir[3], t, x[3], pcoords[3];
  int plane0, plane1, subId;
  vtkMath::Subtract(p2, p1, rayDir);

  if (vtkBox::IntersectWithLine(bounds, p1, p2, t0, t1, x0, x1, plane0, plane1) == 0)
  {
    return 0; // No intersections possible, line is outside the locator
  }

  // Okay process line
  int* ndivs = this->Binner->Divisions;
  vtkIdType prod = ndivs[0] * ndivs[1];
  double* h = this->Binner->H;
  T i, numCellsInBin;
  int ijk[3], ijkEnd[3];
  vtkIdType idx, cId;
  double hitCellBoundsPosition[3], tHitCell;
  double step[3], next[3], tMax[3], tDelta[3];
  double binBounds[6];

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

  // Get the i-j-k point of intersection and bin index. This is
  // clamped to the boundary of the locator. Also get the exit bin
  // of the line.
  this->Binner->GetBinIndices(x0, ijk);
  this->Binner->GetBinIndices(x1, ijkEnd);
  idx = ijk[0] + ijk[1] * ndivs[0] + 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]);