Locate cell containing a point

Development
1

Hi,

This is a duplicate from the VTK discourse. As I didn’t get any answer, I thought it may be more relevant here. Feel free to answer either here or in the VTK discourse, and close the least relevant.

I’m trying to locate a point (with xyz coordinates) in a grid (i.e. find the cell containing that point). It can be any type of grid: image data, structured, explicit structured, unstructured. My goal is to get interpolation weigths to interpolate vertex data.

I’m using the FindCell method of vtkDataSet (https://vtk.org/doc/nightly/html/classvtkDataSet.html#a2221c10d3c4cca44e82c5ef70e4e1cbd , see code below).

It works fine when the dataset is an image data, but for StructuredGrids, I have a crash (see stack trace below).

Am I doing something wrong? Or is there a better way to do what I’m trying to do?

Any help will be greatly appreciated

Thanks

Simon

Code:

// Get the cell containing the point
int sub_id;
double pcoords[3];
vtkIdType cell_id = dataSet->FindCell(point, nullptr, -1, 0, sub_id, pcoords, nullptr);

//Get the number of points of the cell
vtkNew<vtkIdList> cell_points;
dataSet->GetCellPoints(cell_id, cell_points);
vtkIdType nb_points = cell_points->GetNumberOfIds();

//Get the interpolation weights
dataSet->FindCell(point, nullptr, cell_id, 0, sub_id, pcoords, weights);

Stacktrace error:
… my method
7 0x7fe7ebe82c61 vtkPointSet::FindCell(double*, vtkCell*, long long, double, int&, double*, double*) + 99
6 0x7fe7ebe82bb2 vtkPointSet::FindCell(double*, vtkCell*, vtkGenericCell*, long long, double, int&, double*, double*) + 188
5 0x7fe7ebc96afd vtkClosestPointStrategy::FindCell(double*, vtkCell*, vtkGenericCell*, long long, double, int&, double*, double*) + 805
4 0x7fe7ebc967af /home/fuet/SPDE/ParaView/ParaView57/build/lib/libvtkCommonDataModel-pv5.7.so.1(+0x1e97af) [0x7fe7ebc967af]
3 0x7fe7ebc96624 /home/fuet/SPDE/ParaView/ParaView57/build/lib/libvtkCommonDataModel-pv5.7.so.1(+0x1e9624) [0x7fe7ebc96624]
2 0x7fe7ebd13a90 vtkHexahedron::EvaluatePosition(double const*, double*, int&, double*, double&, double*) + 738
1 0x7fe7ebd14382 vtkHexahedron::InterpolationFunctions(double const*, double*) + 114
0 0x7fe7f4153f20 /lib/x86_64-linux-gnu/libc.so.6(+0x3ef20) [0x7fe7f4153f20]
( 29.308s) [paraview ] :0 FATL| Signal: SIGSEGV
Erreur de segmentation (core dumped)

2

Hello,

It seems that the interpolation weights should not be nullptr. I think the following modifications fix the problem.

// Get the cell containing the point
int sub_id;
double pcoords[3];
double weights[8];
//vtkIdType cell_id = dataSet->FindCell(point, nullptr, -1, 0, sub_id, pcoords, nullptr);
vtkIdType cell_id = dataSet->FindCell(point, nullptr, -1, 0, sub_id, pcoords, weights);
3

Hi!

That was simple…
For structured grids it’s ok because every cell has 8 vertices, so weights will always be of size 8.

However for unstructured grids, I don’t know how many points the containing cell has, so I can’t initialize the weights variable, hence my 2 calls to dataSet->FindCell, the first one to retrieve the cell and find out how many points it has, and the second one to compute the interpolation weights.
For now I’m working on structured grids so I’ll stick with this solution. But as in a few weeks I’ll start working with unstructured grids, I’m interested in a more general solution (if there is one).

Thanks a lot!

Simon

4

Hello,

The proper way might be to initialize a large-enough array of weights. Since VTK_CELL_SIZE(=512) is most of the time large enough to represent the number of points contained in an arbitrary cell, this parameter is usually used.

double weights[VTK_CELL_SIZE];

See example:

5

Thanks a lot for your answers! That’s exactly what I need

Simon

6

I’ve added the link to this post in the similar post on the VTK discourse

1 Like