[Insight-users] Bug in itkKdTree
motes motes
mort.motes at gmail.com
Tue May 18 11:09:29 EDT 2010
Here the example with 36 points and without the PointLocator class.
Further I have just updated to itk 3.18.
When this is run only 33 of the 36 points are found even though radius=1000.0:
#ifndef UNIT_TEST_KDTREE_H
#define UNIT_TEST_KDTREE_H
#include "itkPoint.h"
#include "itkPointSet.h"
#include "itkObject.h"
#include "itkKdTree.h"
#include "itkKdTreeGenerator.h"
#ifdef ITK_USE_REVIEW_STATISTICS
#include "itkPointSetToListSampleAdaptor.h"
#else
#include "itkPointSetToListAdaptor.h"
#endif
void unit_test_kdtree() {
typedef float PixelType;
const unsigned int Dimension = 2;
typedef itk::PointSet< PixelType, Dimension >
PointSetType;
typedef PointSetType::PointType PointType;
typedef PointSetType::PointsContainerPointer
PointsContainerPointer;
PointSetType::Pointer PointSet = PointSetType::New();
PointsContainerPointer points = PointSet->GetPoints();
//create 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, p25, p26, p27, p28, p29;
PointType p30, p31, p32, p33, p34, p35;
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]= 50.0; p05[1]= 0.0; points->InsertElement( 5, p05 );
p06[0]= 0.0; p06[1]= 10.0; points->InsertElement( 6, p06 );
p07[0]= 10.0; p07[1]= 10.0; points->InsertElement( 7, p07 );
p08[0]= 20.0; p08[1]= 10.0; points->InsertElement( 8, p08 );
p09[0]= 30.0; p09[1]= 10.0; points->InsertElement( 9, p09 );
p10[0]= 40.0; p10[1]= 10.0; points->InsertElement( 10, p10 );
p11[0]= 50.0; p11[1]= 10.0; points->InsertElement( 11, p11 );
p12[0]= 0.0; p12[1]= 20.0; points->InsertElement( 12, p12 );
p13[0]= 10.0; p13[1]= 20.0; points->InsertElement( 13, p13 );
p14[0]= 20.0; p14[1]= 20.0; points->InsertElement( 14, p14 );
p15[0]= 30.0; p15[1]= 20.0; points->InsertElement( 15, p15 );
p16[0]= 40.0; p16[1]= 20.0; points->InsertElement( 16, p16 );
p17[0]= 50.0; p17[1]= 20.0; points->InsertElement( 17, p17 );
p18[0]= 0.0; p18[1]= 30.0; points->InsertElement( 18, p18 );
p19[0]= 10.0; p19[1]= 30.0; points->InsertElement( 19, p19 );
p20[0]= 20.0; p20[1]= 30.0; points->InsertElement( 20, p20 );
p21[0]= 30.0; p21[1]= 30.0; points->InsertElement( 21, p21 );
p22[0]= 40.0; p22[1]= 30.0; points->InsertElement( 22, p22 );
p23[0]= 50.0; p23[1]= 30.0; points->InsertElement( 23, p23 );
p24[0]= 0.0; p24[1]= 40.0; points->InsertElement( 24, p24 );
p25[0]= 10.0; p26[1]= 40.0; points->InsertElement( 25, p25 );
p26[0]= 20.0; p26[1]= 40.0; points->InsertElement( 26, p26 );
p27[0]= 30.0; p27[1]= 40.0; points->InsertElement( 27, p27 );
p28[0]= 40.0; p28[1]= 40.0; points->InsertElement( 28, p28 );
p29[0]= 50.0; p29[1]= 40.0; points->InsertElement( 29, p29 );
p30[0]= 0.0; p30[1]= 50.0; points->InsertElement( 30, p30 );
p31[0]= 10.0; p31[1]= 50.0; points->InsertElement( 31, p31 );
p32[0]= 20.0; p32[1]= 50.0; points->InsertElement( 32, p32 );
p33[0]= 30.0; p33[1]= 50.0; points->InsertElement( 33, p33 );
p34[0]= 40.0; p34[1]= 50.0; points->InsertElement( 34, p34 );
p35[0]= 50.0; p35[1]= 50.0; points->InsertElement( 35, p35 );
std::cout << "Total points = " << PointSet->GetNumberOfPoints() << std::endl;
/** Type of the PointSet to List Adaptor. */
#ifdef ITK_USE_REVIEW_STATISTICS
typedef itk::Statistics::PointSetToListSampleAdaptor< PointSetType >
SampleAdaptorType;
#else
typedef itk::Statistics::PointSetToListAdaptor< PointSetType >
SampleAdaptorType;
#endif
/** Types fo the KdTreeGenerator */
typedef itk::Statistics::KdTreeGenerator< SampleAdaptorType >
TreeGeneratorType;
typedef typename TreeGeneratorType::Pointer
TreeGeneratorPointer;
typedef typename TreeGeneratorType::KdTreeType TreeType;
typedef typename TreeType::ConstPointer
TreeConstPointer;
typedef typename TreeType::InstanceIdentifierVectorType
InstanceIdentifierVectorType;
TreeConstPointer m_Tree;
typename SampleAdaptorType::Pointer m_SampleAdaptor =
SampleAdaptorType::New();
typename TreeGeneratorType::Pointer m_KdTreeGenerator =
TreeGeneratorType::New();
// Lack of const-correctness in the PointSetAdaptor should be fixed.
m_SampleAdaptor->SetPointSet(PointSet);
m_SampleAdaptor->SetMeasurementVectorSize(Dimension);
m_KdTreeGenerator->SetSample( m_SampleAdaptor );
m_KdTreeGenerator->SetBucketSize( 16 );
m_KdTreeGenerator->Update();
m_Tree = m_KdTreeGenerator->GetOutput();
PointType query;
query[0] = 0.0; query[1] = 0.0;
PointType point( query );
InstanceIdentifierVectorType result;
double radius = 1000.0;
m_Tree->Search( point, radius, result );
int found = result.size();
std::cout << "Points found = " << found << std::endl;
for(unsigned int i = 0; i < result.size(); i++) {
std::cout << "ID = " << result[i] << " pos = " <<
points->GetElement(result[i] ) << std::endl;
}
}
#endif
Am I doing something wrong or is it a bug in the kdTree?
On Tue, May 18, 2010 at 3:24 PM, motes motes <mort.motes at gmail.com> wrote:
> 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?
>
More information about the Insight-users
mailing list