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

[Ricardo Ferrari]

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();
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091228/cbb7fa59/attachment-0001.htm>