Issues locating cells in vtkUnstructuredGrid

Development
1

I have a vtkUnstructuredGrid that can contain 0D, 1D, 2D or 3D cells, and we need to return the cell containing a 3D point given as input. This works fine for 3D cells, but for 0D/1D/2D we need a tolerance because otherwise the point will almost never lie exactly on the element.

We wanted to use one of the available vtkCellLocator classes to accelerate this query but none of them uses the tolerance parameter according to the documentation (the parameter is present on the API, but marked as vtkNotUsed).

So I tried another method, using the FindClosestPointWithinRadius() API, which should support the tolerance parameter properly and also return the cell corresponding to the closest point found. This seems to work in simple cases (e.g. a dataset entirely made of quads that lie on the same plane) but does not return valid results when the dataset has a more complex structure. Even calling FindClosestPoint() with no radius (which should consider the entire dataset) does not return valid results in most cases, which doesn’t make any sense to me.

I found that there’s a FindCell() method on the vtkDataSet class itself (not sure if this one uses an acceleration structure for the location of the cell, but at this point I just wanted it to work) which takes a tolerance and uses it, so I tried that one as well but with very similar results. Points that lie along the edges of the cell seem to be working, but on the rest of the surface it pretty much always fails. The flat dataset always works though.

I want to understand if this a limitation of the VTK library, or if I’m doing something wrong which is causing this to fail.

I have verified that the datasets I create are correct. In the following image you can see both datasets (step 1 and step 15 of the same animation) exported in ParaView. The location of cells works for the “flat” one and not for the other.

datasets
datasets1140×358 240 KB

I would upload the dataset files, but “new users can not upload attachments”.

Thanks in advance,
PaoloG

2

Hello @pgallo, thanks for posting here for the first time.

VTK currently has the following cell locators, vtkCellLocator, vtkStaticCellLocator, vtkCellTreeLocator, vtkModifiedBSPTree, vtkOBBTree,

  1. vtkCellLocator uses an octree as an acceleration structure, and has no multithreaded build step
  2. vtkStaticCellLocator uses a regular grid as an acceleration structure, and has multithreaded build step
  3. vtkCellTreeLocator uses an optimal BSP Tree, and has no multithreaded build step
  4. vtkModifiedBSPTree uses a modified BSP Tree which randomly splits, and has no multithreaded build step
  5. vtkOBBTree is an OBB Tree, and has no multithreaded build step

The first 3 are the most robust out of the 5 cell locators available. vtkCellTreeLocator is the fastest w.r.t. querying, and vtkStaticCellLocator is the fastest w.r.t. build time step. vtkCellLocator and vtkStaticCellLocator are the ones that support the most cell locator related operations.

I provided this information so that you can choose wisely the best locator for you.

Now back to your question.

Tolerance is a vtkLocator parameter, which some vtkAbstractCellLocator/vtkAbstractPointLocator subclasses choose to use or not. For the case of vtkAbstractCellLocator subclasses i have personally included in their header file which parameters do they use or not to improve the documentation. None of the vtkAbstractCellLocator subclasses use the Tolerance parameter, and they are documented accordingly.

vtkAbstractCellLocators have been mostly designed to work for 3D cells

if you would like to use a tolerance parameter for better supporting 0D/1D/2D cells then could you provide a list of the functions that you would like to use? And also provide a dropbox/onedrive/googledrive link with the file that you experimenting and maybe some code too?

Thanks a lot, in advance.

@cory.quammen @will.schroeder

3

Hi @spyridon97, thanks for your reply and the additional explanation of the different locator types. I had actually tried all of them in my attempts at getting this to work so I more or less understood their differences, but I though vtkStaticCellLocator used an octree as well, so I still learned something new.

It makes sense that cell locators are designed to work for 3D cells primarily, and they seem to behave pretty well in that case.

In terms of what functions I need to call, I listed a few in my original post but in practice I’d just need FindCell() to work with a tolerance (either the one from a locator, or just the one from vtkDataSet). The others I mentioned like FindClosestPointWithinRadius() were attempts that I’ve made to work around the issue of locators not supporting a tolerance, but they didn’t work out either.

These are the two datasets that I showed in my example: Step1 (first step of the animation, flat dataset, this works fine), Step15 (a later step of the animation, with non-coplanar cells, this does not work).

Sample source code:

// Find the element which contains the given 3D point.
// In case of 0D/1D/2D elements, it needs a tolerance.
bool VTK_Locator::FindElementFromPoint(vec3f point, double fTolerance, ElementType& eType, sglUInt32 &iElem) const
{
	assert(fTolerance >= 0.0);

#ifdef DISABLE

	if (!_cellLocator)
		return false;

	// NOTE: vtkCellLocator's FindElement API takes a 'tol2' parameter but ignores it.
	// Therefore it is not able to locate the containing element if the point does not
	//  lie EXACTLY inside of it, which is almost impossible for 0D/1D/2D elements in
	//  the floating-point world.
	// To work around this issue, I am instead calling FindClosestPointWithinRadius(),
	//  using the tolerance value as search radius and then returning the element
	//  which contains the closest point (if any).

	double pt_in[3] = { point.x, point.y, point.z };
	double pt_out[3];
	vtkIdType cellId;
	int subId;
	double dist2;
	vtkIdType ret = _cellLocator->FindClosestPointWithinRadius(pt_in, fTolerance, pt_out, cellId, subId, dist2);

#endif//#ifdef DISABLE

	if (!_dataset)
		return false;

	const double fTol2 = fTolerance * fTolerance;
	double pt_in[3] = { point.x, point.y, point.z };
	int subId = 0;
	double weights[8] = { 0.0 };
	double pcoords[3] = { 0.0 };
	vtkIdType cellId = _dataset->FindCell(pt_in, nullptr, -1, fTol2, subId, pcoords, weights);

    /*
      Some code to retrieve element type and ID from our internal data structure, given the VTK cell ID.
    */
    if (!bValidElement)
        return false;
    return true;
}

Thanks for the help,
PaoloG

4

While i look at the implementation details of adding a tolerance for the cell locators for the FindCell operation, could you provide me with a set of sample query points for the step15, that you would expect them to work but they don’t?

5

I have the following VTK MR that you could either check yourself, if it fixes the issue, or you can provide me with points to query for the step15 dataset, and tell you if there is a difference. Be sure to use cellLocator->FindCell function (or the ones lists int the commits of the MR.

https://gitlab.kitware.com/vtk/vtk/-/merge_requests/10578

6

Thanks for working on this @spyridon97.

I am not able to test this right now, so here’s a few sample points that are used for queries at step 15:

{-454.06317138671875, -73.721076965332031, 95.833091735839844}
{-304.56289672851562, -157.73828125000000, 123.16339874267578}
{-342.64373779296875, -3.0045547485351562, 113.86793518066406}
{-540.86804199218750, 19.044498443603516, 83.401062011718750}
{-185.92204284667969, -144.25280761718750, 121.96329498291016}

I think the current tolerance we use is 1.0e-9, but we can go higher. In my tests with FindClosestPointWithinRadius() it was not working even with much larger values like 1.

Regards,
PaoloG

7

I wrote the following code:

    vtkNew<vtkUnstructuredGridReader> reader;
    reader->SetFileName("/home/local/KHQ/spiros.tsalikis/Documents/CellLocators For less than 3d/dataset_step15.vtk");
    reader->Update();
    auto ug = reader->GetOutput();

    std::vector<std::array<double, 3>> samplePoints = {{-454.06317138671875, -73.721076965332031, 95.833091735839844}, {-304.56289672851562, -157.73828125000000, 123.16339874267578}, {-342.64373779296875, -3.0045547485351562, 113.86793518066406}, {-540.86804199218750, 19.044498443603516, 83.401062011718750}, {-185.92204284667969, -144.25280761718750, 121.96329498291016}};
    vtkNew<vtkStaticCellLocator> cellLocator;
    cellLocator->SetDataSet(ug);
    cellLocator->BuildLocator();


    vtkNew<vtkGenericCell> cell;
    double pcoords[3];
    std::vector<double> weights(128, 0);
    vtkIdType cellId;
    int subId;
    double tolerance = 0.0;
    for (auto& samplePoint: samplePoints)
    {
        auto cellId = cellLocator->FindCell(samplePoint.data(), tolerance, cell, subId, pcoords, weights.data());
        std::cout << "cellId: " << cellId << std::endl;
    }

And the results i get with 0 tolerance (so it should be the same as the VTK master branch) are the followings:

cellId: 1946
cellId: 1527
cellId: 2943
cellId: 297
cellId: 1601

Did you get the same results?

Do the results seem correct to you?

8

Currently, it is a bit difficult for me to verify this, for several reasons.

We are currently in “release rush” mode, I couldn’t wait for this to get fixed on the VTK side (also because updating VTK for us is not straight-forward, so it would have taken quite a bit of time). Therefore, I’ve decided to roll my own octree implementation and cell location queries (with tolerance support), which currently seem to be working fine. I’m still using VTK for checking the point against 3D cells (with vtkCell::EvaluatePosition), but I am no longer using any VTK locator in my code.

While implementing this solution, we figured out a significant issue: the query points come from intersecting the geometry with rays, but here the tessellation of a quad cell into 2 triangles was being done differently than “usual” (I consider 0,1,2 / 0,2,3 the “usual” split of a quad into 2 trias, and instead the code was using 0,1,3 / 1,2,3). Since most of these quads are not coplanar, this actually made a difference and it was causing my own location queries to fail in many cases, so I think that a similar issue might have happened with the VTK locator as well (if you’re internally splitting the quad into 2 triangles with 0,1,2 / 0,2,3 indices).

Because of this, those sample points that I gave might not actually be correct with our updated tessellation, and I don’t have the VTK code at hand to test them quickly anyways. I could get back to this once we’re out of this busy period, but unfortunately I’m not able to provide any quick feedback.

With all that said, though, I’m surprised by your results. If I recall correctly, I was not getting any valid cell ID back from the FindCell calls…

9

No worries.

Feel free to get back to me when you have more time, and we will verify that my fixes for tolerance in Cell Locators do what you want them to do.