Understanding vtkStaticCellLocator

Hello Everyone,

I have a case where I need to perform hundreds of thousands or millions IntersectWithLine calls through a cell locator on a geometry dataset that might have millions of cells present in it. In order to perfom it quick. I have made my custom class where I use the vtkCellLocator and vtkSMPTools with TBB backend.

I then tried to compare my code with vtkSelectEnclosedPoints where the function is threaded with SMP and that class uses vtkStaticCellLocator. I have seen comments and posts that vtkStaticCellLocator can improve speed but when I try to use the same. I probably get same result in speed or sometimes worse than the vtkCellLocator. So I want to have an understanding of the Static alternatives of vtk and how they should be used. If there is any relevant documentation. It would be helpful.

Below is an template example for my SMP implementation.

namespace
{
struct SMPClass
{
  float* targets;
  vtkCellLocator* CellLocator;
  vtkPolyData* pointsPolyData;
  
  // SMP Local Object
  vtkSMPThreadLocalObject<vtkGenericCell> genCell;

  SMPClass(float* pointArray,vtkCellLocator* CellLocator,
    vtkPolyData* pointsPolyData)
	: targets(pointArray)
    , CellLocator(CellLocator)
    , pointsPolyData(pointsPolyData)
  {

  }
  ~SMPClass() {}

  void Initialize() {}

  void operator()(vtkIdType startId, vtkIdType endId)
  {
    double point[3];
    vtkGenericCell* genericCell = this->genCell.Local();
	float* target =this->targets + startId;
    for (int i = startId; i < endId; i++)
    {
      this->pointsPolyData->GetPoints()->GetPoint(i, point);
      if (*target) 
      {
        if (CellLocator->IntersectWithLine(point,endPoint,tolerance, nullptr, nullptr, cell))
        {
          // code that modifies target pointer 

        }
      }
	  //Target pointer moves to next element
    }
  }

  void Reduce() {}

  static void Execute(float* pointsptr,vtkCellLocator* CellLocator,
    vtkPolyData* pointsPolyData)
  {
    // Get the instant of the structure
    SMPClass module(pointsptr,CellLocator, pointsPolyData);
    // Use SMP tools for method to execute the operator
    vtkSMPTools::For(0, numPts, module);
  }
};
}


//function that implements Intersectwithline in smp
void customClass::customMethod(vtkFloatArray* pointArray, vtkPolyData* pointsData, vtkCellLocator* cellLocator)
{
  float* pointsptr = static_cast<float*>(pointArray->GetVoidPointer(0)); 
  // Call the static function.
  SMPClass::Execute(pointsptr, cellLocator, pointsData);
}

PS, The cell locator is built outside the class and is only passed to this custom class as a function argument.

Thank you!

1 Like

This does not sound good. By utilizing the massively parallel execution of graphics hardware and/or by using images for computations instead of meshes you may achieve several magnitudes faster computations. Using SMP on CPU you can expect maybe one magnitude faster computations.

Hi Andras,

Thank you for your reply!

I was wondering on how I can implement this using graphical resources or image computation. Could you please give me a recommendation on which classes/libraries to look at. It will be very much helpful.

Thanks once Again!

Have a good day!

You can convert your meshes to an image and use VTK’s GPU volume raycast mapper for processing. You may use custom shaders for computation and the rendered image can return the results. For more abstract computations, you may use compute shaders (for example, you might give this new webGPU compute shader class a try).

Hi Andras,

Thank you for the suggestions! I was able to explore a whole lot about vtk the past few days.

I have been tried your suggestion to convert polydata to image data. However, I have seen that this conversion was very bad as it resulted in some loss of data and I have observed this conversion is slow as well. Here is my code that I have used for conversion:

def vtkPolytoImage(polydata):
    image=vtk.vtkImageData()
    xi, xf, yi, yf, zi, zf = polydata.GetBounds()
    spacing=[0.5,0.5,0.5]
    image.SetSpacing(spacing)
    print(xi, xf, yi, yf, zi, zf)
    #Calculate dimensions
    dx = int(np.ceil((xf - xi) / spacing[0]))
    dy = int(np.ceil((yf - yi) / spacing[1]))
    dz = int(np.ceil((zf - zi) / spacing[2]))
    image.SetDimensions((dx,dy,dz))
    image.SetExtent(0,dx-1,0,dy-1,0,dz-1)
    
    origin=[(x+y)/2 for x,y in zip([xi,yi,zi],spacing)]
    image.SetOrigin(origin)
    image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR,1)
    
    for i in range(image.GetNumberOfPoints()):
        image.GetPointData().GetScalars().SetTuple1(i,255)
    
    pol2stenc=vtk.vtkPolyDataToImageStencil()
    pol2stenc.SetInputData(polydata)
    pol2stenc.SetOutputOrigin(origin)
    pol2stenc.SetOutputSpacing(spacing)
    pol2stenc.SetOutputWholeExtent(image.GetExtent())
    pol2stenc.Update()
    
    imgstenc=vtk.vtkImageStencil()
    imgstenc.SetInputData(image)
    imgstenc.SetStencilConnection(pol2stenc.GetOutputPort())
    imgstenc.ReverseStencilOff()
    imgstenc.SetBackgroundValue(0)
    imgstenc.Update()
    
    writer=vtk.vtkXMLImageDataWriter()
    writer.SetFileName("test.vti")
    writer.SetInputData(imgstenc.GetOutput())
    writer.Write()
    return

I have started learning about shaders after this. To be honest I am completely lost. I don’t understand how this can be used as a replacement of ray casting which I tried to do in cpu earlier. As of my knowlefge, shaders determine how your scene is represented/ rendered onto your screen. Currently my use case is that I want to be able to shoot thousands of rays from a single point to multiple points. I am really stuck on the way forward with this.

Appreciate your help so far!

Thanks!

You can choose a resolution that provides sufficient accuracy for your use case.

This is just a one-time conversion. Once it is done, you may be able to do raycasting magnitudes faster.

If you have a NVIDIA GPU with RT cores then you can try OptiX backend in VTK to perform hardware-accelerated raytracing directly on meshes.

You can have a look at how VTK implements hardware-accelerated raycasting using shaders for volumes in vtkGPUVolumeRayCastMapper and raytracing in Rendering/RayTracing.