Dear Itk-users,<br><br>I am working on the following template function used to estimate initial parameters for a Gaussian Mixture Model. I can compile and link the code without any errors. However, the results I am getting from the filter itk::Statistics::CovarianceSampleFilter are wrong.<br>
<br>Could anybody please check the code to see if I am using this filter properly? <br><br>Thank you,<br>Ricardo<br><br><br><br>/// *****************************************************************************************************<br>
/// <br>/// *****************************************************************************************************<br>template< class TArrayImageType ><br>void TKMeansEstimator< TArrayImageType ><br>::Run()<br>
{<br> ///<br> /// Gather info from the base class<br> ///<br> unsigned int NumberOfClasses = TGMM_InitStrategy< TArrayImageType >::m_NumberOfClasses;<br> unsigned int NumberOfContrasts = TGMM_InitStrategy< TArrayImageType >::m_NumberOfContrasts; <br>
const TArrayImageType *InputImage = TGMM_InitStrategy< TArrayImageType >::m_InputImage;<br><br> ///<br> /// Generate sample data<br> ///<br> typedef itk::Statistics::ImageToListSampleAdaptor< TArrayImageType > ImageToListSampleAdaptorType;<br>
typename ImageToListSampleAdaptorType::Pointer adaptor = ImageToListSampleAdaptorType::New();<br> adaptor->SetImage( InputImage );<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> typedef typename SampleClassifierFilterType::MembershipSampleType::ClassLabelType ClassLabelType;<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>
{<br> ClassLabelType classLabel = (i + 1) * 100;<br> classLabelVector.push_back( classLabel );<br> }<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>
<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> const typename SampleClassifierFilterType::MembershipSampleType* membershipSample = classifier->GetOutput();<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> /// Get a pointer only to the samples of classLabel <br> ClassLabelType classLabel = (i + 1) * 100;<br> const ClassSampleType* classSample = membershipSample->GetClassSample( classLabel );<br>
<br> ///<br> <b> /// Compute covariance matrix <br> /// The problem seems to be here<br> /// <br> typename CovarianceSampleFilterType::Pointer covarianceFilter = CovarianceSampleFilterType::New(); <br>
covarianceFilter->SetInput( classSample );<br> covarianceFilter->Update();<br></b> <br> const typename CovarianceSampleFilterType::MeasurementVectorDecoratedType *MeanDecorator = covarianceFilter->GetMeanOutput();<br>
typename CovarianceSampleFilterType::MeasurementVectorType mean = MeanDecorator->Get();<br> <br> const typename CovarianceSampleFilterType::MatrixDecoratedType *CovDecorator = covarianceFilter->GetCovarianceMatrixOutput();<br>
typename CovarianceSampleFilterType::MatrixType cov = CovDecorator->Get();<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 >( mean[j] ) );<br>
<br> /// Copy covariance matrix<br> for( unsigned k=0; k < NumberOfContrasts; k++ )<br> params->SetElement( index++, static_cast< double >( cov[j][k] ) );<br> } <br>
<br> TGMM_InitStrategy< TArrayImageType >::m_Parameters.push_back( params ); <br> TGMM_InitStrategy< TArrayImageType >::m_WeightClassesVector.SetElement( i, (double)classSample->GetTotalFrequency() / adaptor->GetTotalFrequency() );<br>
<br> }<br>}<br><br> <br>The input image is created as <br><br><br>///<br>const int NumberOfContrasts = 1; //3;<br>const int Dimension = 3;<br><br>/// Image definition<br>typedef signed short InputPixelType;<br>
typedef itk::Image< InputPixelType, Dimension > InputImageType;<br><br>typedef itk::FixedArray< InputPixelType, NumberOfContrasts > ArrayPixelType;<br>typedef itk::Image< ArrayPixelType, Dimension > ArrayImageType;<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>const TArrayImageType *InputImage = caster->GetOutput();<br>
<br><br><br> <br>