[Insight-users] Problems with itk::Statistics::CovarianceSampleFilter

[Karthik Krishnan]

Would you mind posting a minimal compiling test that demonstrates the problem ?

Thanks

On Mon, Dec 28, 2009 at 8:23 PM, Ricardo Ferrari
<rjf.araraquara at gmail.com> wrote:
> Dear Itk-users,
>
> 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.
>
> Could anybody please check the code to see if I am using this filter
> properly?
>
> Thank you,
> Ricardo
>
>
>
> ///
> *****************************************************************************************************
> ///
> ///
> *****************************************************************************************************
> template< class TArrayImageType >
> void TKMeansEstimator< TArrayImageType >
> ::Run()
> {
>     ///
>     /// Gather info from the base class
>     ///
>     unsigned int NumberOfClasses = TGMM_InitStrategy< TArrayImageType
>>::m_NumberOfClasses;
>     unsigned int NumberOfContrasts = TGMM_InitStrategy< TArrayImageType
>>::m_NumberOfContrasts;
>     const TArrayImageType *InputImage = TGMM_InitStrategy< TArrayImageType
>>::m_InputImage;
>
>     ///
>     /// Generate sample data
>     ///
>     typedef itk::Statistics::ImageToListSampleAdaptor< TArrayImageType >
> ImageToListSampleAdaptorType;
>     typename ImageToListSampleAdaptorType::Pointer adaptor =
> ImageToListSampleAdaptorType::New();
>     adaptor->SetImage( InputImage );
>
>     ///
>     /// Create KdTree
>     ///
>     typedef itk::Statistics::WeightedCentroidKdTreeGenerator <
> ImageToListSampleAdaptorType > TreeGeneratorType;
>     typename TreeGeneratorType::Pointer treeGenerator =
> TreeGeneratorType::New();
>     treeGenerator->SetSample( adaptor );
>     treeGenerator->SetBucketSize( m_BuckedSize );
>     treeGenerator->Update();
>
>     typedef itk::Statistics::KdTreeBasedKmeansEstimator < typename
> TreeGeneratorType::KdTreeType > EstimatorType;
>     typename EstimatorType::Pointer estimator = EstimatorType::New();
>
>     ///
>     /// Estimate classes mean vector
>     ///
>     typename EstimatorType::ParametersType initialMeans( NumberOfContrasts *
> NumberOfClasses );
>     estimator->SetParameters( initialMeans );
>     estimator->SetKdTree( treeGenerator->GetOutput() );
>     estimator->SetMaximumIteration( m_MaxNumberOfIterations );
>     estimator->SetCentroidPositionChangesThreshold(
> m_CentroidChangePosThresh );
>     estimator->StartOptimization();
>
>       typename EstimatorType::ParametersType estimatedMeans =
> estimator->GetParameters();
>
>     ///
>     /// Compute classes labels
>     ///
>     typedef itk::Statistics::SampleClassifierFilter<
> ImageToListSampleAdaptorType >  SampleClassifierFilterType;
>     typedef typename
> SampleClassifierFilterType::MembershipSampleType::ClassLabelType
> ClassLabelType;
>
>     typedef typename SampleClassifierFilterType::ClassLabelVectorObjectType
> ClassLabelVectorObjectType;
>     typename ClassLabelVectorObjectType::Pointer classLabelsObject =
> ClassLabelVectorObjectType::New();
>
>     typedef typename SampleClassifierFilterType::ClassLabelVectorType
> ClassLabelVectorType;
>       ClassLabelVectorType& classLabelVector = classLabelsObject->Get();
>     for( unsigned int i=0; i < NumberOfClasses; i++ )
>     {
>         ClassLabelType classLabel = (i + 1) * 100;
>         classLabelVector.push_back( classLabel );
>     }
>
>       /// Set a decision rule type
>       typedef itk::Statistics::MinimumDecisionRule2  DecisionRuleType;
>       typename DecisionRuleType::Pointer decisionRule =
> DecisionRuleType::New();
>
>       const typename
> SampleClassifierFilterType::MembershipFunctionVectorObjectType
> *membershipFunctionsObject = estimator->GetOutput();
>
>       /// Instantiate and pass all the required inputs to the filter
>     typename SampleClassifierFilterType::Pointer classifier =
> SampleClassifierFilterType::New();
>       classifier->SetInput( adaptor );
>       classifier->SetNumberOfClasses( NumberOfClasses );
>       classifier->SetClassLabels( classLabelsObject );
>       classifier->SetDecisionRule( decisionRule );
>       classifier->SetMembershipFunctions( membershipFunctionsObject );
>      classifier->Update();
>     const typename SampleClassifierFilterType::MembershipSampleType*
> membershipSample = classifier->GetOutput();
>
>     ///
>     /// Compute the covariance matrices and the weight vectors
>     ///
>     typedef typename
> SampleClassifierFilterType::MembershipSampleType::ClassSampleType
> ClassSampleType;
>     typedef itk::Statistics::CovarianceSampleFilter< ClassSampleType >
>                 CovarianceSampleFilterType;
>
>     /// Resize weight class vector
>     TGMM_InitStrategy< TArrayImageType >::m_WeightClassesVector.SetSize(
> NumberOfClasses );
>
>     /// Now, get the mean & covariance into the parameters vector
>     for( unsigned int i=0; i < NumberOfClasses; ++i )
>     {
>         /// Get a pointer only to the samples of classLabel
>         ClassLabelType classLabel = (i + 1) * 100;
>         const ClassSampleType* classSample =
> membershipSample->GetClassSample( classLabel );
>
>           ///
>         /// Compute covariance matrix
>         /// The problem seems to be here
>         ///
>         typename CovarianceSampleFilterType::Pointer covarianceFilter =
> CovarianceSampleFilterType::New();
>         covarianceFilter->SetInput( classSample );
>         covarianceFilter->Update();
>
>           const typename
> CovarianceSampleFilterType::MeasurementVectorDecoratedType *MeanDecorator =
> covarianceFilter->GetMeanOutput();
>           typename CovarianceSampleFilterType::MeasurementVectorType mean =
> MeanDecorator->Get();
>
>           const typename CovarianceSampleFilterType::MatrixDecoratedType
> *CovDecorator = covarianceFilter->GetCovarianceMatrixOutput();
>           typename CovarianceSampleFilterType::MatrixType cov =
> CovDecorator->Get();
>
>         /// Create a parameter array for the current class
>         ParametersType *params = new ParametersType( NumberOfContrasts *
> NumberOfContrasts + NumberOfContrasts );
>
>         unsigned int index = 0;
>         for( unsigned j=0; j < NumberOfContrasts; j++ )
>         {
>             /// Copy mean vector
>             params->SetElement( index++, static_cast< double >( mean[j] ) );
>
>             /// Copy covariance matrix
>             for( unsigned k=0; k < NumberOfContrasts; k++ )
>                 params->SetElement( index++, static_cast< double >(
> cov[j][k] ) );
>         }
>
>         TGMM_InitStrategy< TArrayImageType >::m_Parameters.push_back( params
> );
>         TGMM_InitStrategy< TArrayImageType
>>::m_WeightClassesVector.SetElement( i,
> (double)classSample->GetTotalFrequency() / adaptor->GetTotalFrequency() );
>
>     }
> }
>
>
> The input image is created as
>
>
> ///
> const int NumberOfContrasts = 1; //3;
> const int Dimension = 3;
>
> /// Image definition
> typedef signed short
> InputPixelType;
> typedef itk::Image< InputPixelType, Dimension >
> InputImageType;
>
> typedef itk::FixedArray< InputPixelType, NumberOfContrasts >
> ArrayPixelType;
> typedef itk::Image< ArrayPixelType, Dimension >
> ArrayImageType;
>
>
> ///
> /// Create a image of array from different image contrasts
> ///
>     typedef itk::ScalarToArrayCastImageFilter< InputImageType,
> ArrayImageType > CasterType;
>     CasterType::Pointer caster = CasterType::New();
>     for( int i=0; i < NumberOfContrasts; ++i )
>         caster->SetInput( i, ReadMincImage< InputImageType >(
> inputFileName[i] ) );
>     caster->Update();
>
> const TArrayImageType *InputImage = caster->GetOutput();
>
>
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>



-- 
_________________________________
Karthik Krishnan
R&D Engineer,
Kitware Inc.
Ph: +1 5188814919, +91 9538477060