[vtkusers] vtkPointLocator

[Bryn Lloyd]

Hi again,


If anybody is interested in this functionality, I have written the 
functions necessary for finding the N closest points (in point the 
insertion mode), i.e. vtkPointLocator::FindClosestNInsertedPoints and 
vtkPointLocator::FindInsertedPointsWithinRadius.

Maybe these functions could be added to the CVS? Thanks!

Cheers, Bryn




void vtkPointLocator::FindClosestNInsertedPoints(int N, double x, double 
y, double z,
                                   vtkIdList *result)
{
   double xyz[3];
   xyz[0] = x;
   xyz[1] = y;
   xyz[2] = z;
   this->FindClosestNInsertedPoints(N,xyz,result);
}

void vtkPointLocator::FindClosestNInsertedPoints(int N, const double 
x[3], vtkIdList *result)
{
   int i, j;
   double dist2;
   double *pt;
   int level;
   vtkIdType ptId, cno;
   vtkIdList *ptIds;
   int ijk[3], *nei;
   vtkNeighborPoints buckets;

   // clear out the result
   result->Reset();

   /// this->BuildLocator(); // will subdivide if modified; otherwise 
returns

   //
   //  Find bucket point is in.
   //
   for (j=0; j<3; j++)
     {
     ijk[j] = static_cast<int>(
       ((x[j] - this->Bounds[2*j]) /
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

     if (ijk[j] < 0)
       {
       ijk[j] = 0;
       }
     else if (ijk[j] >= this->Divisions[j])
       {
       ijk[j] = this->Divisions[j] - 1;
       }
     }

   // there are two steps, first a simple expanding wave of buckets until
   // we have enough points. Then a refinement to make sure we have the
   // N closest points.
   level = 0;
   double maxDistance = 0.0;
   int currentCount = 0;
   idsort *res = new idsort [N];

   this->GetBucketNeighbors (&buckets, ijk, this->Divisions, level);
   while (buckets.GetNumberOfNeighbors() && currentCount < N)
     {
     for (i=0; i<buckets.GetNumberOfNeighbors(); i++)
       {
       nei = buckets.GetPoint(i);
       cno = nei[0] + nei[1]*this->Divisions[0] +
             nei[2]*this->Divisions[0]*this->Divisions[1];

       if ( (ptIds = this->HashTable[cno]) != NULL )
         {
         for (j=0; j < ptIds->GetNumberOfIds(); j++)
           {
           ptId = ptIds->GetId(j);
           /// pt = this->DataSet->GetPoint(ptId);
           pt = this->Points->GetPoint(ptId);
           dist2 = vtkMath::Distance2BetweenPoints(x,pt);
           if (currentCount < N)
             {
             res[currentCount].dist = dist2;
             res[currentCount].id = ptId;
             if (dist2 > maxDistance)
               {
               maxDistance = dist2;
               }
             currentCount++;
             if (currentCount == N)
               {
               qsort(res, currentCount, sizeof(idsort), vtkidsortcompare);
               }
             }
           else if (dist2 < maxDistance)
             {
             res[N-1].dist = dist2;
             res[N-1].id = ptId;
             qsort(res, N, sizeof(idsort), vtkidsortcompare);
             maxDistance = res[N-1].dist;
             }
           }
         }
       }
     level++;
     this->GetBucketNeighbors (&buckets, ijk, this->Divisions, level);
     }

   // do a sort
   qsort(res, currentCount, sizeof(idsort), vtkidsortcompare);

   // Now do the refinement
   this->GetOverlappingBuckets (&buckets, x, ijk, 
sqrt(maxDistance),level-1);

   for (i=0; i<buckets.GetNumberOfNeighbors(); i++)
     {
     nei = buckets.GetPoint(i);
     cno = nei[0] + nei[1]*this->Divisions[0] +
       nei[2]*this->Divisions[0]*this->Divisions[1];

     if ( (ptIds = this->HashTable[cno]) != NULL )
       {
       for (j=0; j < ptIds->GetNumberOfIds(); j++)
         {
         ptId = ptIds->GetId(j);
         /// pt = this->DataSet->GetPoint(ptId);
         pt = this->Points->GetPoint(ptId);
         dist2 = vtkMath::Distance2BetweenPoints(x,pt);
         if (dist2 < maxDistance)
           {
           res[N-1].dist = dist2;
           res[N-1].id = ptId;
           qsort(res, N, sizeof(idsort), vtkidsortcompare);
           maxDistance = res[N-1].dist;
           }
         }
       }
     }

   // Fill in the IdList
   result->SetNumberOfIds(currentCount);
   for (i = 0; i < currentCount; i++)
     {
     result->SetId(i,res[i].id);
     }

   delete [] res;
}





void vtkPointLocator::FindInsertedPointsWithinRadius(double R, double x, 
double y, double z,
                                       vtkIdList *result)
{
   double xyz[3];
   xyz[0] = x;
   xyz[1] = y;
   xyz[2] = z;
   this->FindInsertedPointsWithinRadius(R,xyz,result);
}

void vtkPointLocator::FindInsertedPointsWithinRadius(double R, const 
double x[3],
                                       vtkIdList *result)
{
   int i, j;
   double dist2;
   double *pt;
   vtkIdType ptId, cno;
   vtkIdList *ptIds;
   int ijk[3], *nei;
   double R2 = R*R;
   vtkNeighborPoints buckets;

   /// this->BuildLocator(); // will subdivide if modified; otherwise 
returns
   //
   //  Find bucket point is in.
   //
   for (j=0; j<3; j++)
     {
     ijk[j] = static_cast<int>(
       ((x[j] - this->Bounds[2*j]) /
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

     if (ijk[j] < 0)
       {
       ijk[j] = 0;
       }
     else if (ijk[j] >= this->Divisions[j])
       {
       ijk[j] = this->Divisions[j] - 1;
       }
     }

   // get all buckets within a distance
   this->GetOverlappingBuckets (&buckets, x, ijk, R, 0);
   // add the original bucket
   buckets.InsertNextPoint(ijk);

   // clear out the result
   result->Reset();

   for (i=0; i<buckets.GetNumberOfNeighbors(); i++)
     {
     nei = buckets.GetPoint(i);
     cno = nei[0] + nei[1]*this->Divisions[0] +
       nei[2]*this->Divisions[0]*this->Divisions[1];

     if ( (ptIds = this->HashTable[cno]) != NULL )
       {
       for (j=0; j < ptIds->GetNumberOfIds(); j++)
         {
         ptId = ptIds->GetId(j);
         /// pt = this->DataSet->GetPoint(ptId);
         pt = this->Points->GetPoint(ptId);
         dist2 = vtkMath::Distance2BetweenPoints(x,pt);
         if (dist2 <= R2)
           {
           result->InsertNextId(ptId);
           }
         }
       }
     }
}









> 
> 
> I would like to use the vtkPointLocator in the point insertion mode, 
> i.e. I would like to incrementally insert points and search for them.
> 
> The incremental point insertion mode works with the functions: 
> InsertNextPoint, IsInsertedPoint, FindClosestInsertedPoint
> 
> 
> There is, however, no method to search for the N closest points.
> 
> I tried FindClosestNPoints with no success. I tried first calling 
> BuildLocator on the locator but get error messages (No points to 
> subdivide).
> What do I need to do, in order to use the function FindClosestNPoints in 
> point insertion mode?
> 
> 
> 
> Another function that is "missing" is point removal, although I can 
> probably live with this.
> 
> Is there a better code for doing what I need?
> 
> 
> 
> Thanks for any suggestions & help.
> 
> Bryn
> 
> _______________________________________________
> This is the private VTK discussion list.
> Please keep messages on-topic. Check the FAQ at: 
> http://www.vtk.org/Wiki/VTK_FAQ
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>