new function: vtkCellLocator::FindCellsWithinRadius(..)
This issue was created automatically from an original Mantis Issue. Further discussion may take place here.
I would like to see a function FindCellsWithinRadius(..) get added to vtkCellLocator (like FindPointsWithinRadius(..) in vtkPointLocator).
Here the source for vtkCellLocator.cxx:
void vtkZFSCellLocator::FindCellsWithinRadius (double R, double x, double y, double z, vtkIdList *result) { double ptr[3]; ptr[0] = x; ptr[1] = y; ptr[2] = z; this->FindCellsWithinRadius(R,ptr,result); }
void vtkZFSCellLocator::FindCellsWithinRadius (double radius, double x[3], vtkIdList *result) { int ijk[3]; vtkIdType j; vtkIdType cellId; vtkIdList *cellIds; double cellBounds[6]; double distance2ToCellBounds;
int leafStart = this->NumberOfOctants- this->NumberOfDivisions*this->NumberOfDivisions*this->NumberOfDivisions;
double radius2 = radius*radius;
double distance2ToDataBounds, maxDistance;
distance2ToDataBounds = this->Distance2ToBounds(x, this->Bounds);
maxDistance = sqrt(distance2ToDataBounds) + this->DataSet->GetLength();
if (radius > maxDistance) {
radius = maxDistance;
radius2 = maxDistance*maxDistance;
}
// Clear the array that indicates whether we have visited this cell.
// The array is only cleared when the query number rolls over. This
// saves a number of calls to memset.
this->QueryNumber++;
if (this->QueryNumber == 0) {
this->ClearCellHasBeenVisited();
this->QueryNumber++; // can't use 0 as a marker
}
// Find bucket point x[3] is in.
//
for (j=0; j<3; j++) {
ijk[j] = (int)((x[j] - this->Bounds[2*j]) / this->H[j]);
if (ijk[j] < 0) {
ijk[j] = 0;
}
else if (ijk[j] >= this->NumberOfDivisions) {
ijk[j] = this->NumberOfDivisions-1;
}
}
// // 1. Start by searching the bucket that the point is in. // if ((cellIds = this->Tree[ leafStart + ijk[0] + ijk[1]*this->NumberOfDivisions + ijk[2]this->NumberOfDivisionsthis->NumberOfDivisions ]) != NULL ) // found cells in this bucket { // query each cell for (j=0; j < cellIds->GetNumberOfIds(); j++) { // get the cell cellId = cellIds->GetId(j); if (this->CellHasBeenVisited[cellId] != this->QueryNumber) { this->CellHasBeenVisited[cellId] = this->QueryNumber;
// check whether we could be close enough to the cell by
// testing the cell bounds
if (this->CacheCellBounds) {
distance2ToCellBounds = this->Distance2ToBounds(x, this->CellBounds[cellId]);
}
else {
this->DataSet->GetCellBounds(cellId, cellBounds);
distance2ToCellBounds = this->Distance2ToBounds(x, cellBounds);
}
if (distance2ToCellBounds < radius2) { // found a cell within radius
result->InsertUniqueId(cellId);
}
} // if (this->CellHasBeenVisited[cellId])
}
}
// // 2.Now, search only those buckets that are within a radius. // // init int i; int *nei;
int prevMinLevel[3], prevMaxLevel[3];
prevMinLevel[0] = prevMaxLevel[0] = ijk[0];
prevMinLevel[1] = prevMaxLevel[1] = ijk[1];
prevMinLevel[2] = prevMaxLevel[2] = ijk[2];
// Build up the list of buckets within the radius
this->GetOverlappingBuckets(x, ijk, radius, prevMinLevel, prevMaxLevel);
for (i=0; i<this->Buckets->GetNumberOfNeighbors(); i++) {
nei = this->Buckets->GetPoint(i);
if ( (cellIds = this->Tree[leafStart
+ nei[0]
+ nei[1]*this->NumberOfDivisions
+ nei[2]*this->NumberOfDivisions*this->NumberOfDivisions //nei[2]*numberOfBucketsPerPlane
]) != NULL ) {
for (j=0; j < cellIds->GetNumberOfIds(); j++) {
// get the cell
cellId = cellIds->GetId(j);
if (this->CellHasBeenVisited[cellId] != this->QueryNumber) {
this->CellHasBeenVisited[cellId] = this->QueryNumber;
// check whether we could be close enough to the cell by
// testing the cell bounds
if (this->CacheCellBounds) {
distance2ToCellBounds = this->Distance2ToBounds(x, this->CellBounds[cellId]);
}
else {
this->DataSet->GetCellBounds(cellId, cellBounds);
distance2ToCellBounds = this->Distance2ToBounds(x, cellBounds);
}
if (distance2ToCellBounds < radius2) { // found a cell within radius
result->InsertUniqueId(cellId);
}
} // if (this->CellHasBeenVisited[cellId])
}
}
}
}