[vtkusers] vtkPointLocator::FindPointsWithinRadius NOT thread safe ?

[Jan Ehrhardt]

Hi Folks,

we have some trouble using multithreading and
vtkPointLocator::FindPointsWithinRadius(), VTK version 5.4.2.

All is fine, if we run FindPointsWithinRadius() from a single thread.
But running our code in parallel finds incorrect points.
BuildLocator() is called before multithreading.

I include some source code. The code works well with numberOfThreads =
1. If numberOfThreads > 1 lots of Error-Messages "Error: Distance >
Radius: ..." appear.

Thanks for your help,

Jan 

void TestMultithreadding( vtkPolyData *surface )
{
  mihINFO("  Build point locators ...");

  unsigned int numPoints = surface->GetNumberOfPoints();

  // build the locator
  vtkSmartPointer<vtkPointLocator> spPointLocator =
vtkSmartPointer<vtkPointLocator>::New();
  spPointLocator->SetDataSet( surface );
  spPointLocator->Initialize();
  spPointLocator->BuildLocator();

  const unsigned int numberOfThreads = 4;
  int ret;
  int thr_error = 0;
  unsigned int thread_loop;
  pthread_t process_id[24];
  pthread_attr_t attr;

  pthread_attr_init( &attr );

#ifndef __sgi
  pthread_attr_setscope( &attr, PTHREAD_SCOPE_SYSTEM );
#endif

  mihINFO("  Start Threads ("<<numberOfThreads<<")...");

  // Test if for all threads a function is set.
  for( thread_loop = 0; thread_loop < numberOfThreads; thread_loop++ )
  {
    // ==========================================
    // Calculate smoothed points.
    mihSmoothPointsParameters* forceParameters;
    unsigned int threadRange = numPoints / numberOfThreads;
    for( unsigned int threadNr = 0; threadNr < numberOfThreads;
++threadNr )
    {
      // Initialize thread parameters.
      forceParameters = new mihSmoothPointsParameters();

      // Calculate thread range.
      forceParameters->from = threadNr * threadRange;
      forceParameters->to = (threadNr == numberOfThreads - 1) ?
numPoints : (threadNr + 1) * threadRange;

      // Initializing thread parameters.
      forceParameters->surface = surface;
      forceParameters->pointLocator = spPointLocator;
      forceParameters->numPoints = numPoints;
      forceParameters->maxNumOfThreads = numberOfThreads;
      forceParameters->threadNum = thread_loop;

      if( thread_loop < numberOfThreads - 1 )
      {
        ret = pthread_create( &(process_id[thread_loop]), &attr,
T_SmoothPoints, (void *) forceParameters );

        if( ret != 0 )
        {
          mihERROR( "Execution of pthread_create failed !" );
        }
      }
      // Start the last thread-method.
      else
      {
        T_SmoothPoints( (void *) forceParameters );
      }
    }
  }

  // The parent thread has finished its method - so now it waits for
each of the
  // other processes (created with sproc) to exit.
  for( thread_loop = 0; thread_loop < numberOfThreads - 1; thread_loop
++ )
  {
    int ret = pthread_join( process_id[thread_loop], NULL );
    if( ret != 0 )
    {
      mihERROR( "Error waiting for the termination of a thread
(pthread_join)!" );
      thr_error = THR_WAIT_ERROR;
    }
  }
  mihINFO("  Finished.");

}

///////////////////////////////////////////////////////////////////////////////////////////
//
// T_SmoothPoints()
//
///////////////////////////////////////////////////////////////////////////////////////////
MIH_THREAD_RETURN_TYPE T_SmoothPoints( void* vargs )
{
  mihSmoothPointsParameters* args = (mihSmoothPointsParameters*) vargs;

  vtkSmartPointer<vtkPolyData> surface = args->surface;
  vtkSmartPointer<vtkPointLocator> pointLocator = args->pointLocator;

  const double radius = 5;

  double surfacePoint1[3];
  double surfacePoint2[3];

  double dist2 = 0.0;

  // Define neighborhood ID list
  vtkIdList *idList = vtkIdList::New();

  //std::cout << "From " << args->from << " To " << args->to <<
std::endl;
  for( unsigned int i = args->from; i < args->to; ++i )
  {
    surface->GetPoint( i, surfacePoint1 );

    idList->Reset();
    pointLocator->FindPointsWithinRadius( radius, surfacePoint1,
idList );

    int foundPoints = idList->GetNumberOfIds();
    if( foundPoints < 1 )
    {
      mihERROR( "No Points in Radius " << radius << "found !" );
      continue;
    }

    for( int j = 0; j < foundPoints; j++ )
    {
      int id = idList->GetId( j );

      if( id >= static_cast<int> ( args->numPoints ) || id < 0 )
      {
        mihERROR( "Invalid Id: " << id );
        continue;
      }

      surface->GetPoint( id, surfacePoint2 );

      dist2 = vtkMath::Distance2BetweenPoints( surfacePoint1,
surfacePoint2 );
      if( dist2 > radius2 )
      {
        mihERROR( "Distance > Radius: " << dist2 );
        continue;
      }
    }
  }

  // Delete the list
  idList->Delete();

  return MIH_THREAD_RETURN_VALUE;
}