vtkStaticCellLocator problem (vtk9.0.3)

marcomusy · February 16, 2022, 5:40pm

I observe something strange in vtkStaticCellLocator.FindClosestPoint() vs vtkCellLocator.FindClosestPoint():

with vtkStaticCellLocator (wrong):

with vtkCellLocator (correct):

from vedo import *

# Load and cut mouse limb mesh
limb = Mesh("https://vedo.embl.es/examples/data/270_flank.vtk")
limb.alpha(0.3).c("purple9").lighting("off")

# Make a wall in the cut limb to show gradient
grid = Grid(sx=2000, sy=2000, resx=50, resy=50).rotateX(90)
grid.cutWithMesh(limb).wireframe(False)

# scaling things makes no difference
#limb.scale(0.001)
#grid.scale(0.001)

values = []
for p in grid.points():
    q = limb.closestPoint(p)
    values.append(mag(p-q))
grid.cmap("viridis_r", values).addScalarBar()

show(limb, grid, viewup='z', axes=1)

The range is also very different.
The name of the two classes is literally the only thing I change between the two… Any ideas?

Yohann_Bearzi · February 16, 2022, 5:56pm

This seems to be a bug. What happens if you run the same script but replace vtkStaticCellLocator by vtkStaticPointLocator?

will.schroeder · February 16, 2022, 6:01pm

The behavior of the two classes should be the same. Could you create an equivalent VTK python script…

marcomusy · February 16, 2022, 6:55pm

…unless i’m doing something silly…:

import vtk
from vedo import Mesh, Grid, mag, show

# load the data as VTK objs
limb = Mesh("https://vedo.embl.es/examples/data/270_flank.vtk")
lpoints = limb.points()  # numpy array
limb = limb.polydata()   # vtkPolyData
grid = Grid(sx=2000, sy=2000, resx=50, resy=50)
gpoints = grid.points()  # numpy array
grid = grid.polydata()   # vtkPolyData

locator1 = vtk.vtkCellLocator()
locator1.SetDataSet(limb)
locator1.BuildLocator()

locator2 = vtk.vtkStaticCellLocator()
locator2.SetDataSet(limb)
locator2.BuildLocator()

locator3 = vtk.vtkPointLocator()
locator3.SetDataSet(limb)
locator3.BuildLocator()

locator4 = vtk.vtkStaticPointLocator()
locator4.SetDataSet(limb)
locator4.BuildLocator()

values1 = []
values2 = []
values3 = []
values4 = []
for p in gpoints:
    q = [0, 0, 0]
    cid = vtk.mutable(0)
    dist2 = vtk.mutable(0)
    subid = vtk.mutable(0)
    locator1.FindClosestPoint(p, q, cid, subid, dist2)
    values1.append(mag(p-q))

    q = [0, 0, 0]
    cid = vtk.mutable(0)
    dist2 = vtk.mutable(0)
    subid = vtk.mutable(0)
    locator2.FindClosestPoint(p, q, cid, subid, dist2)
    values2.append(mag(p-q))

    id3 = locator3.FindClosestPoint(p)
    values3.append(mag(p-lpoints[id3]))

    id4 = locator4.FindClosestPoint(p)
    values4.append(mag(p-lpoints[id3]))

print(max(values1), max(values2)) # should be the same
print(max(values3), max(values4)) # should be the same (and they are)

# show(Mesh(grid).cmap("jet", values2))

outputs:

795.8835103198636 993.0246205222754
795.8835103198636 795.8835103198636

Screenshot from 2022-02-16 19-39-49

Yohann_Bearzi · February 16, 2022, 7:07pm

I just run this script in ParaView master and I do get the correct max values. What version of VTK are you using?

marcomusy · February 16, 2022, 7:09pm

it’s
vtk version : 9.0.3
python version : 3.8.8 (default, Apr 13 2021, 19:58:26) [GCC 7.3.0]
python interpreter: /home/musy/soft/anaconda3/bin/python
system : Linux 5.4.0-99-generic posix x86_64

will.schroeder · February 16, 2022, 7:12pm

There was a bug fix made to vtkStaticCellLocator related to FindClosestPoint (commit 2b07d64ea made about a year ago…)

Yohann_Bearzi · February 16, 2022, 7:21pm

Yep, you should update to VTK 9.1.0 to have this fixed.

Jens_Munk_Hansen · March 16, 2023, 3:50pm

I stumbled upon the same issue, when replacing the vtkCellLocator of vtkImplicitPolyDataDistance with the vtkStaticCellLocator to get some extra performance. Is there a reason why the vtkStaticCellLocator is not used as the default in vtkImplicitPolyDataDistance. My experience is that neither vtkCellLocator or vtkStaticCellLocator are thread-safe. Would it be an idea to introduce thread_local into the hierarchy of class for vtkStaticCellLocator to make it thread-safe? Making it thread-safe would make it possible to use vtkSMPTools are bit more widely, e.g. for the vtkPolyDataDistanceFilter?

will.schroeder · March 17, 2023, 11:09am

Is there a reason why the vtkStaticCellLocator is not used as the default in vtkImplicitPolyDataDistance

Probably because vtkStaticCellLocator was added after vtkImplicitPolyDataDistance, and vtkImplicitPolyDataDistance was never updated.

vtkStaticCellLocator should be thread safe (once built). There is an outstanding issue for a potential race condition that needs tracking down, maybe this is biting you?

Jens_Munk_Hansen · March 17, 2023, 4:28pm

Thanks. I get some weird indices but I am also not waiting for the locators to build. I build a separate locator for each thread in the Initialize function using vtkSMPTools. But I guess from the name, vtkStaticCellLocator, this obviously goes wrong.

I will wait for a single locator to build and use it across multiple threads. I am sure this should work.

If I track down and locate a data race, I will ensure to communicate it such that errors are fixed. Thanks for replying

Jens

will.schroeder · March 17, 2023, 5:18pm

Generally what I’ve done is build the locator prior to the parallel for. Mostly because the locators can chew up a lot of memory, and can take some time to build, so I minimize these by using just the one instance.

The race is reported by Sean who is (rightfully) bugging me

Jens_Munk_Hansen · March 17, 2023, 5:20pm

Exactly what I am doing now and things are looking pretty fast here - more memory efficient than what I did before.

What I did before, which most likely triggers the data race was the following.

Replace vtkCellLocator with vtkStaticCellLocator in vtkImplicitPolyDataDistance and use the functor below.

class MyOp
{
private:
  vtkSMPThreadLocalObject<vtkImplicitPolyDataDistance> Imp;
  vtkSpsDisplacementPolyDataFilter* Self;
  vtkDoubleArray* OutputScalar;
  vtkDoubleArray* OutputVector;
  vtkPolyData* Src;
  vtkPolyData* Dst;
public:
  MyOp(vtkDisplacementPolyDataFilter* self,
       vtkImplicitPolyDataDistance* src, vtkPolyData* dst,
       vtkDoubleArray* outputScalar,
       vtkDoubleArray* outputVector) :
    Self(self), Src(src), Dst(dst),
    OutputScalar(outputScalar), OutputVector(outputVector)
  {
    
  }
  void Initialize()
  {
    // Weird things happen
    vtkSpsImplicitPolyDataDistance*& imp = this->Imp.Local();
    imp->SetInput(this->Src);
  }
  void operator()(vtkIdType begin, vtkIdType end)
  {
    vtkImplicitPolyDataDistance*& imp = this->Imp.Local();
    for (vtkIdType ptId = begin; ptId < end; ptId++) {
      double pt[3];
      this->Dst->GetPoint(ptId, pt);
    
      double closestPoint[3];
      double val = imp->EvaluateFunctionAndGetClosestPoint(pt, closestPoint);
      double dist = Self->SignedDistance ? (Self->NegateDistance ? -val : val) : fabs(val);
      double direction[3];
      double absDistance = 0.0;
      for (int i = 0 ; i < 3 ; i++) {
        direction[i] = closestPoint[i] - pt[i];
      }
      absDistance = std::max<double>(1e-6, fabs(dist));

      // Normalize direction
      for (int i = 0 ; i < 3 ; i++) {
        direction[i] /= absDistance;
      }
      
      OutputVector->SetTuple(ptId, direction);
      OutputScalar->SetValue(ptId, dist);
    }
  }
  void Reduce() {
  }
};

I know that I can use the gradient which is faster than computing a new vector and normalizing it… Thanks again for reaching out.

levi · May 13, 2023, 9:03am

Hi, Jens

i have just replace vtkCellLocator with vtkStaticCellLocator in vtkImplicitPolyDataDistance, with nothing else.
but the process is as slow as vtkCellLocator,
i`m sure smp is supported since other classes with smp is quite fast
am i missing something? what other change should i take?
thanks

Jens_Munk_Hansen · May 13, 2023, 10:22am

I have replaced vtkCellLocator and introduced vtkSMP. It made it 2 orders of magnitude faster. Below is a snippet from a new vtkDisplacementFilter, where both distance and directions are computed in parallell.

namespace
{
class DisplacementOp
{
private:
  vtkDisplacementPolyDataFilter* Self;
  vtkDoubleArray* OutputScalar;
  vtkDoubleArray* OutputVector;
  vtkImplicitPolyDataDistance* Imp;
  vtkPolyData* Dst;

public:
  DisplacementOp(vtkDisplacementPolyDataFilter* self, vtkImplicitPolyDataDistance* imp,
                 vtkPolyData* dst, vtkDoubleArray* outputScalar, vtkDoubleArray* outputVector)
    : Self(self)
    , Imp(imp)
    , Dst(dst)
    , OutputScalar(outputScalar)
    , OutputVector(outputVector) {}
  void Initialize() {}
  void operator()(vtkIdType begin, vtkIdType end)
  {
    vtkImplicitPolyDataDistance* imp = this->Imp;
    vtkTypeBool signedDistance = Self->GetSignedDistance();
    vtkTypeBool negateDistance = Self->GetNegateDistance();
    for (vtkIdType ptId = begin; ptId < end; ptId++)
      {
        double pt[3];
        Dst->GetPoint(ptId, pt);

        double closestPoint[3];
        double val = imp->EvaluateFunctionAndGetClosestPoint(pt, closestPoint);
        double dist = signedDistance
           ? (negateDistance ? -val : val) : fabs(val);
        double direction[3];
        double absDistance = 0.0;
        for (int i = 0; i < 3; i++)
          {
            direction[i] = closestPoint[i] - pt[i];
          }
        absDistance = std::max<double>(1e-6, fabs(dist));

        // Normalize direction
        for (int i = 0; i < 3; i++)
          {
            direction[i] /= absDistance;
          }

        OutputVector->SetTuple(ptId, direction);
        OutputScalar->SetValue(ptId, dist);
      }
  }
  void Reduce() {}
};

and to execute it

  vtkImplicitPolyDataDistance* imp = vtkImplicitPolyDataDistance::New();
  imp->SetInput(src);

  // Calculate distance from points.
  int numPts = mesh->GetNumberOfPoints();

  vtkDoubleArray* directionArray = vtkDoubleArray::New();
  directionArray->SetName("Direction");
  directionArray->SetNumberOfComponents(3);
  directionArray->SetNumberOfTuples(numPts);

  vtkDoubleArray* distanceArray = vtkDoubleArray::New();
  distanceArray->SetName("Distance");
  distanceArray->SetNumberOfComponents(1);
  distanceArray->SetNumberOfTuples(numPts);

  // The sign of the functions are the dot product between pseudo normal and the vector
  // x-p. The point p is on the target polydata.
  DisplacementOp pointFunctor(this, imp, mesh, distanceArray, directionArray);
  vtkSMPTools::For(0, numPts, pointFunctor);

levi · May 15, 2023, 9:05am

Hi Jens
thanks for your reply, i have changed vtkDistancePolyDataFilter like yours and it works.
though the time is just reduced about 75%(i will continue optimize it).

with vtk 9.1, i have some error result in vtkStaticCellLoacator, i’ll have some check

by the way, i found your commit in vtk for https://gitlab.kitware.com/vtk/vtk/-/commit/90cf6e9cd3876bf7423af7f4ed00df3ce620ef21

Jens_Munk_Hansen · May 15, 2023, 9:53pm

I got a factor of around 70 but in addition to the commit you mentioned, I replaced the vtkCellLocator with vtkStaticCellLocator. Also remember to enable vtkSMPTools. I usually use TBB.

levi · May 16, 2023, 1:16am

got it
by now i just use OpenMP, maybe i should test against TBB

Jens_Munk_Hansen · May 16, 2023, 3:05pm

or even just STDThread.