<br>Dear ITK-Users,<br><br>I am working on a class to estimate the initial
parameters (mean vector, covariance matrix and weighting vector) for a
Gaussian Mixture model. For that I am using the
KdTreeBasedKmeansEstimator and the refactoring statistic classes. To be
more specifically, I am using the itk::Statistics::<div id=":4n" class="ii gt">CovarianceSampleFilter
to compute the Mean vector and Covariance matrices. However, the
results I am getting are completely wrong. After tracking the code, it
seems that the mean vector is always zero. <br>
<br>Could anybody help me to solve this problem please? I am not sure if I am missing something. <br><br>Thank you,<br>Ricardo<br><br><br> ///<br> /// Create a image of array from different image contrasts<br> ///<br>
typedef itk::ScalarToArrayCastImageFilter< InputImageType, ArrayImageType > CasterType;<br> CasterType::Pointer caster = CasterType::New();<br> for( int i=0; i < NumberOfContrasts; ++i )<br> caster->SetInput( i, ReadMincImage< InputImageType >( inputFileName[i] ) );<br>
caster->Update();<br><br><br> ///<br> /// Generate sample data<br> ///<br> typedef itk::Statistics::ImageToListSampleAdaptor< TArrayImageType > ImageToListSampleAdaptorType;<br> typename ImageToListSampleAdaptorType::Pointer adaptor = ImageToListSampleAdaptorType::New();<br>
adaptor->SetImage( caster->GetOutput() );<br><br> ///<br> /// Create KdTree<br> /// <br> typedef itk::Statistics::WeightedCentroidKdTreeGenerator < ImageToListSampleAdaptorType > TreeGeneratorType;<br>
typename TreeGeneratorType::Pointer treeGenerator = TreeGeneratorType::New();<br> treeGenerator->SetSample( adaptor );<br> treeGenerator->SetBucketSize( m_BuckedSize );<br> treeGenerator->Update();<br>
<br> typedef itk::Statistics::KdTreeBasedKmeansEstimator < typename TreeGeneratorType::KdTreeType > EstimatorType;<br> typename EstimatorType::Pointer estimator = EstimatorType::New();<br><br> ///<br> /// Estimate classes mean vector<br>
///<br> typename EstimatorType::ParametersType initialMeans( NumberOfContrasts * NumberOfClasses );<br> estimator->SetParameters( initialMeans );<br> estimator->SetKdTree( treeGenerator->GetOutput() );<br>
estimator->SetMaximumIteration( m_MaxNumberOfIterations );<br> estimator->SetCentroidPositionChangesThreshold( m_CentroidChangePosThresh );<br> estimator->StartOptimization();<br><br> typename EstimatorType::ParametersType estimatedMeans = estimator->GetParameters();<br>
<br> ///<br> /// Compute classes labels <br> /// <br> typedef itk::Statistics::SampleClassifierFilter< ImageToListSampleAdaptorType > SampleClassifierFilterType;<br><br> typedef typename SampleClassifierFilterType::ClassLabelVectorObjectType ClassLabelVectorObjectType;<br>
typename ClassLabelVectorObjectType::Pointer classLabelsObject = ClassLabelVectorObjectType::New();<br> <br> typedef typename SampleClassifierFilterType::ClassLabelVectorType ClassLabelVectorType;<br> ClassLabelVectorType& classLabelVector = classLabelsObject->Get();<br>
for( unsigned int i=0; i < NumberOfClasses; i++ )<br> classLabelVector.push_back( typename SampleClassifierFilterType::ClassLabelType( (i + 1) * 100 ) );<br><br> //Set a decision rule type<br> typedef itk::Statistics::MinimumDecisionRule2 DecisionRuleType;<br>
typename DecisionRuleType::Pointer decisionRule = DecisionRuleType::New();<br><br> const typename SampleClassifierFilterType::MembershipFunctionVectorObjectType *membershipFunctionsObject = estimator->GetOutput();<br>
const typename SampleClassifierFilterType::MembershipFunctionVectorType membershipFunctions = membershipFunctionsObject->Get();<br><br> //Instantiate and pass all the required inputs to the filter<br> typename SampleClassifierFilterType::Pointer classifier = SampleClassifierFilterType::New();<br>
classifier->SetInput( adaptor );<br> classifier->SetNumberOfClasses( NumberOfClasses );<br> classifier->SetClassLabels( classLabelsObject );<br> classifier->SetDecisionRule( decisionRule );<br> classifier->SetMembershipFunctions( membershipFunctionsObject );<br>
classifier->Update();<br><br> const typename SampleClassifierFilterType::MembershipSampleType* membershipSample = classifier->GetOutput();<br><br>/*<br> const typename SampleClassifierFilterType::MembershipSampleType::ClassLabelType label1( 100 );<br>
const typename SampleClassifierFilterType::MembershipSampleType::ClassSampleType* subSample1 = membershipSample->GetClassSample( label1 );<br> std::cout << subSample1->GetTotalFrequency() << std::endl;<br>
<br> const typename SampleClassifierFilterType::MembershipSampleType::ClassLabelType label2( 200 );<br> const typename SampleClassifierFilterType::MembershipSampleType::ClassSampleType* subSample2 = membershipSample->GetClassSample( label2 );<br>
std::cout << subSample2->GetTotalFrequency() << std::endl;<br><br> const typename SampleClassifierFilterType::MembershipSampleType::ClassLabelType label3( 300 );<br> const typename SampleClassifierFilterType::MembershipSampleType::ClassSampleType* subSample3 = membershipSample->GetClassSample( label3 );<br>
std::cout << subSample3->GetTotalFrequency() << std::endl;<br>*/ <br> <br> <br><br> ///<br> /// Compute the covariance matrices and the weight vectors<br> ///<br> typedef typename SampleClassifierFilterType::MembershipSampleType::ClassSampleType ClassSampleType;<br>
typedef itk::Statistics::CovarianceSampleFilter< ClassSampleType > CovarianceSampleFilterType;<br> <br> /// Resize weight class vector<br> TGMM_InitStrategy< TArrayImageType >::m_WeightClassesVector.SetSize( NumberOfClasses );<br>
<br> /// Now, get the mean & covariance into the parameters vector<br> for( unsigned int i=0; i < NumberOfClasses; ++i )<br> {<br> typename SampleClassifierFilterType::MembershipSampleType::ClassLabelType classLabel( (i + 1) * 100 );<br>
<br> /// Get a pointer only to the samples of classLabel <br> const ClassSampleType* sampleClass = membershipSample->GetClassSample( classLabel );<br><br>/* The result here is okay<br>std::cout << sampleClass->Size() << std::endl; <br>
<br> double sum = 0;<br> for( int ii=0; ii < sampleClass->Size(); ii++ )<br> {<br> typename ClassSampleType::MeasurementVectorType val = sampleClass->GetMeasurementVectorByIndex( ii );<br>
sum = sum + val[0];<br> }<br>std::cout << "Mean value = " << (double)( sum / sampleClass->Size() ) << std::endl;<br>*/<br><br><br> /// Compute covariance matrix <br>
typename CovarianceSampleFilterType::Pointer covarianceFilter = CovarianceSampleFilterType::New(); <br> covarianceFilter->SetInput( sampleClass );<br> covarianceFilter->Update();<br> <br>
const typename CovarianceSampleFilterType::MeasurementVectorDecoratedType *MeanDecorator = covarianceFilter->GetMeanOutput();<br> typename CovarianceSampleFilterType::MeasurementVectorType mean = MeanDecorator->Get();<br>
<br>/* The results here are wrong */ <br>std::cout << "Mean Vector1: " << MeanDecorator->Get() << std::endl;<br>std::cout << "Mean Vector2: " << covarianceFilter->GetMean() << std::endl;<br>
<br> const typename CovarianceSampleFilterType::MatrixDecoratedType *CovDecorator = covarianceFilter->GetCovarianceMatrixOutput();<br> typename CovarianceSampleFilterType::MatrixType cov = CovDecorator->Get();<br>
<br>/* The results here are wrong */ <br>
std::cout << "Covariance Matrix1: " << CovDecorator->Get() << std::endl;<br>std::cout << "Covariance Matrix2: " << covarianceFilter->GetCovarianceMatrix() << std::endl;<br>
<br> /// Create a parameter array for the current class <br> ParametersType *params = new ParametersType( NumberOfContrasts * NumberOfContrasts + NumberOfContrasts );<br><br> unsigned int index = 0;<br>
for( unsigned j=0; j < NumberOfContrasts; j++ )<br> {<br> /// Copy mean vector<br> params->SetElement( index++, static_cast< double >( estimatedMeans[i * NumberOfContrasts + j] ) );<br>
<br> /// Copy covariance matrix<br> for( unsigned k=0; k < NumberOfContrasts; k++ )<br> params->SetElement( index++, (double)cov[j][k] );<br> } <br><br>
TGMM_InitStrategy< TArrayImageType >::m_Parameters.push_back( params ); <br> TGMM_InitStrategy< TArrayImageType >::m_WeightClassesVector.SetElement( i, (double)sampleClass->GetTotalFrequency() / adaptor->GetTotalFrequency() );<br>
<br> }</div>