[Insight-users] A question about Spectral clustering using KdTreeBasedKMeansClustering.cxx

[johnjtchen at aol.co.uk]

Dear All,

I am using KdTreeBasedKMeansClustering.cxx from example directory to do 
the spectral clustering.
The problem I had was that no matter what input weight matrix I fed in, 
the centroid vectors obtained from 
KdTreeBasedKmeansEstimator::GetParameters() have only only vector (in 
fact the first centroid vector) valid.
Other centroid vectors are filled with 0 and thus the clustering result 
is not correct.
So what steps I did wrong?


I made following changes to the KdTreeBasedKMeansClustering.cxx.

The input image is eigen vector image 
(eigenvec01h26m40s25020801h26m40s250208.nii.nii.gz).
Eigen vectors were stored in columns and the first eigen vector was 
uninformative so we discard it.
The result of clustering should be able to separate the eigen vectors 
into 2 classes.
Each class should contains half-ish eigen vectors.

//////////////////////////////////////////////////////////////////////////
   int numberofclass = 2;

   typedef itk::Image <float,3> ImageType;
   typedef itk::ImageFileReader <ImageType> ReaderType;
   ReaderType::Pointer reader = ReaderType::New();
   reader->SetFileName(inFileName);
   reader->Update();

   ImageType::Pointer weightMatrix = reader->GetOutput();
//////////////////////////////////////////////////////////////////////////

   //     typedef itk::Vector< double, 1 > MeasurementVectorType;
   //     typedef itk::Statistics::ListSample< MeasurementVectorType > 
SampleType;
   //     SampleType::Pointer sample = SampleType::New();
   //     sample->SetMeasurementVectorSize( 1 );
   int measurementVecLen = 
reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
   typedef itk::Array <float> MeasurementVectorType;
   MeasurementVectorType mv;
   mv.SetSize(measurementVecLen);
//////////////////////////////////////////////////////////////////////////

   //     typedef itk::Statistics::NormalVariateGenerator 
NormalGeneratorType;
   //     NormalGeneratorType::Pointer normalGenerator = 
NormalGeneratorType::New();
   //
   //     normalGenerator->Initialize( 101 );
   //
   //     MeasurementVectorType mv;
   //     double mean = 100;
   //     double standardDeviation = 30;
   //     for ( unsigned int i = 0 ; i < 100 ; ++i )
   //     {
   //         mv[0] = ( normalGenerator->GetVariate() * 
standardDeviation ) + mean;
   //         sample->PushBack( mv );
   //     }
   //
   //     normalGenerator->Initialize( 3024 );
   //     mean = 200;
   //     standardDeviation = 30;
   //     for ( unsigned int i = 0 ; i < 100 ; ++i )
   //     {
   //         mv[0] = ( normalGenerator->GetVariate() * 
standardDeviation ) + mean;
   //         sample->PushBack( mv );
   //     }
   typedef itk::Statistics::ListSample <MeasurementVectorType> SampleType;
   SampleType::Pointer sample = SampleType::New();
   sample->SetMeasurementVectorSize(mv.GetSize());
   for ( unsigned int j = 1 ; j < measurementVecLen ; ++j )
   {
       mv.Fill(0);
       for ( unsigned int i = 0 ; i < measurementVecLen ; ++i )
       {
           ImageType::IndexType index;
           index[0] = i;
           index[1] = j;
           index[2] = 0;
           mv[i] = weightMatrix->GetPixel(index);

       }
       std::cout << " mv : " << mv << std::endl;
       sample->PushBack(mv);
   }
//////////////////////////////////////////////////////////////////////////

   //     EstimatorType::ParametersType initialMeans(2);
   //     initialMeans[0] = 0.0;
   //     initialMeans[1] = 0.0;
   EstimatorType::ParametersType initialMeans(measurementVecLen * 
numberofclass);
   initialMeans.fill(0);
//////////////////////////////////////////////////////////////////////////

   //     for ( unsigned int i = 0 ; i < 2 ; ++i )
   //     {
   //         std::cout << "cluster[" << i << "] " << std::endl;
   //         std::cout << "    estimated mean : " << estimatedMeans[i] 
<< std::endl;
   //     }
   int indexformeanvector = 0;
   for ( unsigned int i = 0 ; i < numberofclass ; ++i )
   {
       std::cout << "cluster[" << i << "] " << std::endl;
       for ( unsigned int k = 0 ; k < sample->GetMeasurementVectorSize() 
; ++k ){
           std::cout << "    estimated mean : " << 
estimatedMeans[indexformeanvector++] << std::endl;
       }
   }
//////////////////////////////////////////////////////////////////////////

   std::vector< unsigned int > classLabels;
   //     classLabels.resize( 2 );
   //     classLabels[0] = 100;
   //     classLabels[1] = 200;
   classLabels.resize(numberofclass);
   for (int i=1;i<=classLabels.size();i++)
   {
       classLabels[i-1] = i*100;
   }

//////////////////////////////////////////////////////////////////////////

The full code and test images can be found in

http://rapidshare.com/files/94784236/questionaboutkmeansclustering.ZIP.html


Thanks very much.

john