Examples/Statistics/WeightedSampleStatistics.cxx

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

// Software Guide : BeginLatex
//
// \index{itk::Statistics::Weighted\-Mean\-Calculator}
// \index{itk::Statistics::Weighted\-Covariance\-Calculator}
// \index{Statistics!Weighted mean}
// \index{Statistics!Weighted covariance}
//
// We include the header file for the \doxygen{Vector} class that will
// be our measurement vector template in this example.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkVector.h"
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// We will use the \subdoxygen{Statistics}{ListSample} as our sample
// template. We include the header for the class too.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkListSample.h"
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
//
// The following headers are for the weighted covariance algorithms.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkWeightedMeanSampleFilter.h"
#include "itkWeightedCovarianceSampleFilter.h"
// Software Guide : EndCodeSnippet

using MeasurementVectorType = itk::Vector< float, 3 >;

class ExampleWeightFunction :
  public itk::FunctionBase< MeasurementVectorType, double >
{
public:
  using Self = ExampleWeightFunction;
  using Superclass = itk::FunctionBase< MeasurementVectorType, double >;
  using Pointer = itk::SmartPointer<Self>;
  using ConstPointer = itk::SmartPointer<const Self>;

  itkTypeMacro(ExampleWeightFunction, FunctionBase);
  itkNewMacro(Self);

  using InputType = MeasurementVectorType;

  using OutputType = double;

  OutputType Evaluate( const InputType& input ) const override
    {
      if ( input[0] < 3.0 )
        {
        return 0.5;
        }
      else
        {
        return 0.01;
        }
    }

protected:
  ExampleWeightFunction() = default;
  ~ExampleWeightFunction() override = default;
}; // end of class

int main()
{
  // Software Guide : BeginLatex
  //
  // The following code snippet will create a ListSample instance
  // with three-component float measurement vectors and put five
  // measurement vectors in the ListSample object.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using SampleType = itk::Statistics::ListSample< MeasurementVectorType >;
  SampleType::Pointer sample = SampleType::New();
  sample->SetMeasurementVectorSize( 3 );
  MeasurementVectorType mv;
  mv[0] = 1.0;
  mv[1] = 2.0;
  mv[2] = 4.0;

  sample->PushBack( mv );

  mv[0] = 2.0;
  mv[1] = 4.0;
  mv[2] = 5.0;
  sample->PushBack( mv );

  mv[0] = 3.0;
  mv[1] = 8.0;
  mv[2] = 6.0;
  sample->PushBack( mv );

  mv[0] = 2.0;
  mv[1] = 7.0;
  mv[2] = 4.0;
  sample->PushBack( mv );

  mv[0] = 3.0;
  mv[1] = 2.0;
  mv[2] = 7.0;
  sample->PushBack( mv );
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // Robust versions of covariance algorithms require
  // weight values for measurement vectors. We have two ways of
  // providing weight values for the weighted mean and weighted
  // covariance algorithms.
  //
  // The first method is to plug in an array of weight values. The
  // size of the weight value array should be equal to that of the
  // measurement vectors. In both algorithms, we use the
  // \code{SetWeights(weights)}.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using WeightedMeanAlgorithmType =
    itk::Statistics::WeightedMeanSampleFilter<SampleType>;

  WeightedMeanAlgorithmType::WeightArrayType weightArray( sample->Size() );
  weightArray.Fill( 0.5 );
  weightArray[2] = 0.01;
  weightArray[4] = 0.01;

  WeightedMeanAlgorithmType::Pointer weightedMeanAlgorithm =
                                              WeightedMeanAlgorithmType::New();

  weightedMeanAlgorithm->SetInput( sample );
  weightedMeanAlgorithm->SetWeights( weightArray );
  weightedMeanAlgorithm->Update();

  std::cout << "Sample weighted mean = "
            << weightedMeanAlgorithm->GetMean() << std::endl;

  using WeightedCovarianceAlgorithmType =
    itk::Statistics::WeightedCovarianceSampleFilter<SampleType>;

  WeightedCovarianceAlgorithmType::Pointer weightedCovarianceAlgorithm =
                                        WeightedCovarianceAlgorithmType::New();

  weightedCovarianceAlgorithm->SetInput( sample );
  weightedCovarianceAlgorithm->SetWeights( weightArray );
  weightedCovarianceAlgorithm->Update();

  std::cout << "Sample weighted covariance = " << std::endl;
  std::cout << weightedCovarianceAlgorithm->GetCovarianceMatrix() << std::endl;
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // The second method for computing weighted statistics is to plug-in a
  // function that returns a weight value that is usually a function of each
  // measurement vector. Since the \code{weightedMeanAlgorithm} and
  // \code{weightedCovarianceAlgorithm} already have the input sample plugged
  // in, we only need to call the \code{SetWeightingFunction(weights)}
  // method.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  ExampleWeightFunction::Pointer weightFunction = ExampleWeightFunction::New();

  weightedMeanAlgorithm->SetWeightingFunction( weightFunction );
  weightedMeanAlgorithm->Update();

  std::cout << "Sample weighted mean = "
            << weightedMeanAlgorithm->GetMean() << std::endl;

  weightedCovarianceAlgorithm->SetWeightingFunction( weightFunction );
  weightedCovarianceAlgorithm->Update();

  std::cout << "Sample weighted covariance = " << std::endl;
  std::cout << weightedCovarianceAlgorithm->GetCovarianceMatrix();

  std::cout << "Sample weighted mean (from WeightedCovarainceSampleFilter) = "
            << std::endl << weightedCovarianceAlgorithm->GetMean()
            << std::endl;
  // Software Guide : EndCodeSnippet

  return EXIT_SUCCESS;
}