[Insight-users] Bug in itkKdTree

[motes motes]

I use the itkKdtree through the PointLocator2 class found at:

http://svn.na-mic.org/NAMICSandBox/trunk/QuadEdgeMeshRigidRegistration/Source/

I am using InsightToolkit-3.16.0 build for release on Ubuntu 9.04.


I have manually created 25 2D points ranging from [0;40]. I then use
the point 0.0 as my query point and call the Search function
specifying a radius of 1000.0.

This should give me all points in the container but for some reason
point10 (0, 20) is not returned and only 24 points are found. I have
tried to extend the example to 6*6=36 points and then I only get 33
returned even though I specify a radius of 1000.0.

Below is the code for 25 points:


  typedef float PixelType;
  const unsigned int Dimension = 2;
  typedef itk::PointSet< PixelType, Dimension >
PointSetType;
  typedef PointSetType::PointType                                     PointType;
  typedef PointSetType::PointsContainerPointer
PointsContainerPointer;
  typedef itk::PointLocator2<PointSetType>
PointLocatorType;
  typedef PointLocatorType::InstanceIdentifierVectorType
InstanceIdentifierVectorType;

  PointSetType::Pointer  PointSet = PointSetType::New();
  PointsContainerPointer  points = PointSet->GetPoints();
  PointLocatorType::Pointer locator = PointLocatorType::New();

  // create 25 points
  PointType p00, p01, p02, p03, p04, p05, p06, p07, p08, p09;
  PointType p10, p11, p12, p13, p14, p15, p16, p17, p18, p19;
  PointType p20, p21, p22, p23, p24;

  p00[0]=   0.0; p00[1]= 0.0;     points->InsertElement( 0,  p00 );
  p01[0]=  10.0; p01[1]= 0.0;     points->InsertElement( 1,  p01 );
  p02[0]=  20.0; p02[1]= 0.0;     points->InsertElement( 2,  p02 );
  p03[0]=  30.0; p03[1]= 0.0;     points->InsertElement( 3,  p03 );
  p04[0]=  40.0; p04[1]= 0.0;     points->InsertElement( 4,  p04 );

  p05[0]=   0.0; p05[1]= 10.0;    points->InsertElement( 5,  p05 );
  p06[0]=  10.0; p06[1]= 10.0;    points->InsertElement( 6,  p06 );
  p07[0]=  20.0; p07[1]= 10.0;    points->InsertElement( 7,  p07 );
  p08[0]=  30.0; p08[1]= 10.0;    points->InsertElement( 8,  p08 );
  p09[0]=  40.0; p09[1]= 10.0;    points->InsertElement( 9,  p09 );

  p10[0]=   0.0; p10[1]= 20.0;    points->InsertElement( 10, p10 );
  p11[0]=  10.0; p11[1]= 20.0;    points->InsertElement( 11, p11 );
  p12[0]=  20.0; p12[1]= 20.0;    points->InsertElement( 12, p12 );
  p13[0]=  30.0; p13[1]= 20.0;    points->InsertElement( 13, p13 );
  p14[0]=  40.0; p14[1]= 20.0;    points->InsertElement( 14, p14 );

  p15[0]=   0.0; p15[1]= 30.0;    points->InsertElement( 15, p15 );
  p16[0]=  10.0; p16[1]= 30.0;    points->InsertElement( 16, p16 );
  p17[0]=  20.0; p17[1]= 30.0;    points->InsertElement( 17, p17 );
  p18[0]=  30.0; p18[1]= 30.0;    points->InsertElement( 18, p18 );
  p19[0]=  40.0; p19[1]= 30.0;    points->InsertElement( 19, p19 );

  p20[0]=   0.0; p20[1]= 40.0;    points->InsertElement( 20, p20 );
  p21[0]=  10.0; p21[1]= 40.0;    points->InsertElement( 21, p21 );
  p22[0]=  20.0; p22[1]= 40.0;    points->InsertElement( 22, p22 );
  p23[0]=  30.0; p23[1]= 40.0;    points->InsertElement( 23, p23 );
  p24[0]=  40.0; p24[1]= 40.0;    points->InsertElement( 24, p24 );

  locator->SetPointSet( PointSet );
  locator->Initialize();
  std::cout << "PointSet->GetNumberOfPoints() : " <<
PointSet->GetNumberOfPoints() << std::endl;

  PointType query2;
  query2[0] = 0.0; query2[1] = 0.0;
  PointLocatorType::PointType point2( query2 );
  InstanceIdentifierVectorType result2;
  locator->Search( point2, 1000.0, result2 );

  int found = result2.size();
  std::cout << "points found = " << found << std::endl;
  for(unsigned int i = 0; i < result2.size(); i++) {
    std::cout << "ID = " << result2[i] << " pos = " <<
points->GetElement(result2[i] ) << std::endl;
  }

In the above loop all points are printed except point with ID=10. Any
ideas for a fix for this or should I try to found a non-itk kd-tree?