[vtkusers] Determine which points are inside of a cell

[David Gobbi]

On Mon, Apr 12, 2010 at 10:10 AM, David Doria <daviddoria+vtk at gmail.com> wrote:
> I have a vtkImageData (named 'grid') superimposed on a vtkPolyData
> (named 'polydata').
>
> I can determine which cell a particular point is in like this:
>
>  vtkSmartPointer<vtkCellLocator> cellLocator =
>      vtkSmartPointer<vtkCellLocator>::New();
>  cellLocator->SetDataSet(grid);
>  cellLocator->BuildLocator();
>
>  //determine which cell each point is in
>  cout << "Which cell does each point belong to?" << endl;
>  for(vtkIdType i = 0; i < polydata->GetNumberOfPoints(); i++)
>    {
>    double p[3];
>    polydata->GetPoint(i,p);
>    int cellId = cellLocator->FindCell(p);
>    cout << "Point " << i << " is in cell " << cellId << endl;
>    }
>
> Now I'd like to do the reverse - determine which points are in a
> particular cell. Is there a good way to do this? Or would I have to
> determine which cell every point is in and construct a table, then use
> that table for this second type of query? I could also make a new
> polydata from the single cell and then use vtkSelectEnclosedPoints,
> but that also seems awkward.

Hi David,

Don't use locators with structured data like vtkImageData.  The
FindPoint() and FindCell() methods for the vtkImageData object itself
work in constant time and are hundreds of times faster than using
FindPoint/Cell on a locator.  Also, building a locator for an image
data is wasteful... the locator will eat up around 24 times as much
memory as the image itself (assuming a 16-bit image) because the
locator computes and stores the bounding box for each and every cell.

As an aside, you also have to be careful about what vtk means by a
"vtkVoxel" cell.  The vtkVoxel is a cell with a data point at each
corner, and is used by VTK's generic data interface for interpolating
between the data values at the eight corner points.  It's actually the
points in the vtkImageData that correspond to what people usually mean
when they say "voxel".  As an illustration of this, if you display an
image with vtkImageActor, then the center each square pixel you see on
the screen corresponds to the (x,y,z) of a point in the data.

Another way of explaining this is that when you store data in the
image, you store the data as "PointData", not as "CellData".  So
please make sure that image->FindCell is what you want, because I
suspect that image->FindPoint is the method that actually gives you
the information you need.


Finally getting to your question, about the only really useful method
for finding multiple points is this:

vtkPointLocator::FindPointsWithinRadius(radius, point, pointIdList);

It finds all the points within a sphere, so you need "radius" to be
larger than the voxel so that it will return more points than you
need, and then you can use a "for" loop to exclude points that aren't
within the voxel.  It's clunky, but I don't think there is any other
way.

  David