[Insight-users] Submission: PCA decomposition calculator
Zachary Pincus
zpincus at stanford.edu
Mon, 3 May 2004 13:42:33 -0700
--Apple-Mail-2-83054363
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
OK, here is a new version of the code and the test for the PCA
decomposition calculator.
The test was slightly broken (my fault) which lead to all of the
errors, but now that should be fixed. I've also made the testing a bit
more complete.
I've also tightened the code up a bit now, for what that's worth.
Thanks, and I am quite impressed at how well the kitware build/test
system works.
Zach
--Apple-Mail-2-83054363
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
x-mac-type=54455854;
x-unix-mode=0755;
name="itkImagePCADecompositionCalculator.h"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculator.h
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkImagePCADecompositionCalculator.h,v $
Language: C++
Date: $Date: 2004/05/03 01:33:36 $
Version: $Revision: 1.2 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __itkImagePCADecompositionCalculator_h
#define __itkImagePCADecompositionCalculator_h
#include "itkObject.h"
#include "itkImagePCAShapeModelEstimator.h"
#include "vnl/vnl_vector.h"
#include "vnl/vnl_matrix.h"
namespace itk
{
/** \class ImagePCADecompositionCalculator
* \brief Decomposes an image into directions along basis components.
*
* This calculator computes the projection of an image into a subspace specified
* by some orthonormal basis.
* Typically, this basis will be the principal components of an image data set,
* as calculated by an ImagePCAShapeModelEstimator. The output of the calculator
* is a vnl_vector containing the coefficients along each dimension of the
* provided basis set.
* To use this calculator, first set each basis image with the SetBasisImage
* method. In the PCA case, the basis images are the outputs of the
* ImagePCAShapeModelEstimator (except the zeroth output, which is the average
* image).
* SetBasisFromModel is a convenience method to set all of this information from
* a given ImagePCAShapeModelEstimator instance.
*
* This class is templated over the input image type and the type of images
* used to describe the basis.
*
* \warning This method assumes that the input image consists of scalar pixel
* types.
*
* \warning All images (input and basis) must be the same size.
*
* \author Zachary Pincus
*
* \ingroup Operators
*/
template <class TInputImage,
class TBasisImage = Image<double, ::itk::GetImageDimension<TInputImage>::ImageDimension> >
class ITK_EXPORT ImagePCADecompositionCalculator : public Object
{
public:
/** Standard class typedefs. */
typedef ImagePCADecompositionCalculator Self;
typedef Object Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(ImagePCADecompositionCalculator, Object);
/** Type definitions for the input images. */
typedef TInputImage InputImageType;
typedef TBasisImage BasisImageType;
/** Pointer types for the image. */
typedef typename TInputImage::Pointer InputImagePointer;
typedef typename TBasisImage::Pointer BasisImagePointer;
/** Const Pointer type for the image. */
typedef typename TInputImage::ConstPointer InputImageConstPointer;
/** Basis image pixel type: this is also the type of the optput vector */
typedef typename TBasisImage::PixelType BasisPixelType;
/** Input Image dimension */
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
/** Basis Image dimension */
itkStaticConstMacro(BasisImageDimension, unsigned int,
TBasisImage::ImageDimension);
/** Vector of basis image pointers. */
typedef std::vector< BasisImagePointer > BasisImagePointerVector;
/** Type definitions for internal vectors and matrices */
typedef vnl_matrix<BasisPixelType> BasisMatrixType;
typedef vnl_vector<BasisPixelType> BasisVectorType;
/** Set the input image. */
itkSetConstObjectMacro(Image,InputImageType);
/** Set the basis images. */
void SetBasisImages(const BasisImagePointerVector _arg);
/** Get the basis images. */
BasisImagePointerVector& GetBasisImages(void) {return m_BasisImages;}
/** Type definition of a compatible ImagePCAShapeModelEstimator */
typedef typename ImagePCAShapeModelEstimator<TInputImage,
TBasisImage>::Pointer ModelPointerType;
/** Set the basis images from a ImagePCAShapeModelEstimator */
void SetBasisFromModel(ModelPointerType model);
/** Compute the PCA decomposition of the input image. */
void Compute(void);
/** Return the projection of the image. */
itkGetMacro(Projection,BasisVectorType);
protected:
ImagePCADecompositionCalculator();
virtual ~ImagePCADecompositionCalculator() {};
void PrintSelf(std::ostream& os, Indent indent) const;
void CalculateBasisMatrix(void);
void CalculateImageAsVector(void);
private:
typedef typename BasisImageType::SizeType BasisSizeType;
ImagePCADecompositionCalculator(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
BasisVectorType m_Projection;
BasisVectorType m_ImageAsVector;
BasisImagePointerVector m_BasisImages;
BasisSizeType m_Size;
InputImageConstPointer m_Image;
BasisMatrixType m_BasisMatrix;
bool m_BasisMatrixCalculated;
unsigned long m_NumPixels;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkImagePCADecompositionCalculator.txx"
#endif
#endif
--Apple-Mail-2-83054363
Content-Transfer-Encoding: 7bit
Content-Type: application/text;
x-mac-type=54455854;
x-unix-mode=0755;
name="itkImagePCADecompositionCalculator.txx"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculator.txx
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkImagePCADecompositionCalculator.txx,v $
Language: C++
Date: $Date: 2004/05/03 14:28:07 $
Version: $Revision: 1.2 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itkImagePCADecompositionCalculator_txx
#define _itkImagePCADecompositionCalculator_txx
#include "itkImagePCADecompositionCalculator.h"
#include "itkImageRegionConstIterator.h"
namespace itk
{
/*
* Constructor
*/
template<class TInputImage, class TBasisImage>
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::ImagePCADecompositionCalculator()
{
m_Image = NULL;
m_BasisMatrixCalculated = false;
}
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::SetBasisImages(const BasisImagePointerVector _arg)
{
//itkDebugMacro("setting BasisImages to " << _arg); // Doesn't seem to work!
if (m_BasisImages.size() != _arg.size() ||
!std::equal(m_BasisImages.begin(), m_BasisImages.end(), _arg.begin()))
{
this->m_BasisImages = _arg;
this->Modified();
this->m_BasisMatrixCalculated = false;
// We need this modified setter function so that the calculator
// can cache the basis set between calculations. Note that computing the
// basis matrix from the input images is rather expensive, and the basis
// images are likely to be changed less often than the input images. So
// it makes sense to try to cache the pre-computed matrix.
}
}
/*
* Compute the projection
*/
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::Compute(void)
{
if (!m_BasisMatrixCalculated)
{
this->CalculateBasisMatrix();
}
this->CalculateImageAsVector();
m_Projection = m_BasisMatrix * m_ImageAsVector;
}
/*
* Convert a vector of basis images into a matrix. Each image is flattened into 1-D.
*/
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::CalculateBasisMatrix(void) {
m_Size = m_BasisImages[0]->GetRequestedRegion().GetSize();
m_NumPixels = 1;
for( unsigned int i = 0; i < BasisImageDimension; i++ )
{
m_NumPixels *= m_Size[i];
}
m_BasisMatrix = BasisMatrixType(m_BasisImages.size(), m_NumPixels);
int i = 0;
for(typename BasisImagePointerVector::const_iterator basis_it = m_BasisImages.begin();
basis_it != m_BasisImages.end(); ++basis_it)
{
if( (*basis_it)->GetRequestedRegion().GetSize() != m_Size)
{
itkExceptionMacro("All basis images must be the same size!");
}
ImageRegionConstIterator<BasisImageType> image_it(*basis_it,
(*basis_it)->GetRequestedRegion());
int j = 0;
for (image_it.GoToBegin(); !image_it.IsAtEnd(); ++image_it)
{
m_BasisMatrix(i, j++) = image_it.Get();
}
i++;
}
m_BasisMatrixCalculated = true;
m_ImageAsVector.set_size(m_NumPixels);
}
/*
* Convert an image into a 1-D vector, changing the pixel type if necessary.
*/
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::CalculateImageAsVector(void) {
if ( m_Image->GetRequestedRegion().GetSize() != m_Size)
{
itkExceptionMacro("Input image must be the same size as the basis images!");
}
ImageRegionConstIterator<InputImageType> image_it(m_Image,
m_Image->GetRequestedRegion());
typename BasisVectorType::iterator vector_it;
for (image_it.GoToBegin(), vector_it = m_ImageAsVector.begin();
!image_it.IsAtEnd(); ++image_it, ++vector_it)
{
*vector_it = static_cast<BasisPixelType> (image_it.Get());
}
}
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::SetBasisFromModel(ModelPointerType model) {
BasisImagePointerVector images;
model->Update(); // Is this a good idea?
unsigned int nImages = model->GetNumberOfPrincipalComponentsRequired();
images.reserve(nImages);
for(int i = 1; i <= nImages; i++)
{
images.push_back(model->GetOutput(i));
}
this->SetBasisImages(images);
}
template<class TInputImage, class TBasisImage>
void
ImagePCADecompositionCalculator<TInputImage, TBasisImage>
::PrintSelf( std::ostream& os, Indent indent ) const
{
Superclass::PrintSelf(os,indent);
os << indent << "Projection: " << m_Projection << std::endl;
os << indent << "Image: " << m_Image.GetPointer() << std::endl;
}
} // end namespace itk
#endif
--Apple-Mail-2-83054363
Content-Transfer-Encoding: 7bit
Content-Type: application/text;
x-mac-type=54455854;
x-unix-mode=0755;
name="itkImagePCADecompositionCalculatorTest.cxx"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculatorTest.cxx
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkImagePCADecompositionCalculatorTest.cxx,v $
Language: C++
Date: $Date: 2004/05/02 19:27:54 $
Version: $Revision: 1.1 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
// Insight classes
#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkLightProcessObject.h"
#include "itkTextOutput.h"
#include "itkImagePCADecompositionCalculator.h"
// class to support progress feeback
class ShowProgressObject
{
public:
ShowProgressObject(itk::LightProcessObject * o)
{m_Process = o;}
void ShowProgress()
{std::cout << "Progress " << m_Process->GetProgress() << std::endl;}
itk::LightProcessObject::Pointer m_Process;
};
int main(int, char* [] )
{
itk::OutputWindow::SetInstance(itk::TextOutput::New().GetPointer());
//Data definitions
const unsigned int IMGWIDTH = 2;
const unsigned int IMGHEIGHT = 2;
const unsigned int NDIMENSION = 2;
const unsigned int NUMTRAINIMAGES = 3;
const unsigned int NUMLARGESTPC = 2;
//------------------------------------------------------
//Create 3 simple test images with
//------------------------------------------------------
typedef itk::Image<double,NDIMENSION> InputImageType;
typedef
itk::ImageRegionIterator< InputImageType > InputImageIterator;
InputImageType::Pointer image1 = InputImageType::New();
InputImageType::Pointer image2 = InputImageType::New();
InputImageType::Pointer image3 = InputImageType::New();
InputImageType::Pointer image4 = InputImageType::New();
InputImageType::Pointer image5 = InputImageType::New();
InputImageType::Pointer image6 = InputImageType::New();
InputImageType::Pointer image7 = InputImageType::New();
InputImageType::SizeType inputImageSize = {{ IMGWIDTH, IMGHEIGHT }};
InputImageType::IndexType index;
index.Fill(0);
InputImageType::RegionType region;
region.SetSize( inputImageSize );
region.SetIndex( index );
//--------------------------------------------------------------------------
// Set up Image 1 first
//--------------------------------------------------------------------------
image1->SetRegions( region );
image1->Allocate();
// setup the iterators
InputImageIterator image1It( image1, image1->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 2 first
//--------------------------------------------------------------------------
image2->SetRegions( region );
image2->Allocate();
// setup the iterators
InputImageIterator image2It( image2, image2->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 3 first
//--------------------------------------------------------------------------
image3->SetRegions( region );
image3->Allocate();
// setup the iterators
InputImageIterator image3It( image3, image3->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 4 first
//--------------------------------------------------------------------------
image4->SetRegions( region );
image4->Allocate();
// setup the iterators
InputImageIterator image4It( image4, image4->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 5 first
//--------------------------------------------------------------------------
image5->SetRegions( region );
image5->Allocate();
// setup the iterators
InputImageIterator image5It( image5, image5->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 6 first
//--------------------------------------------------------------------------
image6->SetRegions( region );
image6->Allocate();
// setup the iterators
InputImageIterator image6It( image6, image6->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 7 first
//--------------------------------------------------------------------------
image7->SetRegions( region );
image7->Allocate();
// setup the iterators
InputImageIterator image7It( image7, image7->GetBufferedRegion() );
//--------------------------------------------------------------------------
//Manually create and store each vector
//--------------------------------------------------------------------------
// The first two vectors are the first two principal components of the data:
// [1 1 1 1] , [2 0 0 2], [0 3 3 0]
// The second two vectors are some of those data, which we will project down
// to the PC-space
// The last three vectors are a new basis set to projet down into, so we can
// test changing bases mid-stream.
//Image no. 1
image1It.Set( -0.3853 ); ++image1It;
image1It.Set( 0.5929 ); ++image1It;
image1It.Set( 0.5929 ); ++image1It;
image1It.Set( -0.3853 ); ++image1It;
//Image no. 2
image2It.Set( -0.5929 ); ++image2It;
image2It.Set( -0.3853 ); ++image2It;
image2It.Set( -0.3853 ); ++image2It;
image2It.Set( -0.5929 ); ++image2It;
//Image no. 3
image3It.Set( 2 ); ++image3It;
image3It.Set( 0 ); ++image3It;
image3It.Set( 0 ); ++image3It;
image3It.Set( 2 ); ++image3It;
//Image no. 4
image4It.Set( 0 ); ++image4It;
image4It.Set( 3 ); ++image4It;
image4It.Set( 3 ); ++image4It;
image4It.Set( 0 ); ++image4It;
//Image no. 5
image5It.Set( 0.70710678 ); ++image5It; // 1/sqrt(2)
image5It.Set( 0.70710678 ); ++image5It;
image5It.Set( 0 ); ++image5It;
image5It.Set( 0 ); ++image5It;
//Image no. 6
image6It.Set( -0.70710678 ); ++image6It;
image6It.Set( 0.70710678 ); ++image6It;
image6It.Set( 0 ); ++image6It;
image6It.Set( 0 ); ++image6It;
//Image no. 7
image7It.Set( 0 ); ++image7It;
image7It.Set( 0 ); ++image7It;
image7It.Set( 1 ); ++image7It;
image7It.Set( 0 ); ++image7It;
//----------------------------------------------------------------------
// Test code for the Decomposition Calculator
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//Set the image Decomposition Calculator
//----------------------------------------------------------------------
typedef itk::ImagePCADecompositionCalculator<InputImageType>
ImagePCAShapeModelEstimatorType;
ImagePCAShapeModelEstimatorType::Pointer
decomposer = ImagePCAShapeModelEstimatorType::New();
//----------------------------------------------------------------------
//Set the parameters of the clusterer
//----------------------------------------------------------------------
// add the first two vectors to the projection basis
ImagePCAShapeModelEstimatorType::BasisImagePointerVector basis;
basis.push_back(image1);
basis.push_back(image2);
decomposer->SetBasisImages(basis);
// compute some projections!
ImagePCAShapeModelEstimatorType::BasisVectorType proj3, proj4;
decomposer->SetImage(image3);
decomposer->Compute();
proj3 = decomposer->GetProjection();
decomposer->SetImage(image4);
decomposer->Compute();
proj4 = decomposer->GetProjection();
// get the basis images
ImagePCAShapeModelEstimatorType::BasisImagePointerVector basis_check;
basis_check = decomposer->GetBasisImages();
basis.clear();
basis.push_back(image5);
basis.push_back(image6);
basis.push_back(image7);
decomposer->SetBasisImages(basis);
ImagePCAShapeModelEstimatorType::BasisVectorType proj3_2, proj4_2;
//decomposer->SetImage(image4); // DON'T set image4 -- it should still
// be cached with the decomposer. Test that this works between basis changes.
decomposer->Compute();
proj4_2 = decomposer->GetProjection();
decomposer->SetImage(image3);
decomposer->Compute();
proj3_2 = decomposer->GetProjection();
// get the basis images
ImagePCAShapeModelEstimatorType::BasisImagePointerVector basis_check_2;
basis_check_2 = decomposer->GetBasisImages();
//Test the printself function to increase coverage
decomposer->Print(std::cout);
// Print the basis and projections: first the PCA basis
std::cout << "The basis of projection is: " << std::endl;
for (ImagePCAShapeModelEstimatorType::BasisImagePointerVector::const_iterator
basis_it = basis_check.begin(); basis_it != basis_check.end(); ++basis_it)
{
std::cout << "[";
InputImageIterator basisImage_it( *basis_it, (*basis_it)->GetBufferedRegion() );
for (basisImage_it.GoToBegin(); !basisImage_it.IsAtEnd(); ++basisImage_it)
{
std::cout << basisImage_it.Get() << " ";
}
std::cout << "]" << std::endl;
}
//Print the projections
std::cout << "The projection of [0 2 2 0] is [" << proj3 << "]" << std::endl;
std::cout << "this should be approx [-1.5412 -2.3716]" << std::endl;
std::cout << "The projection of [0 3 3 0] is [" << proj4 << "]" << std::endl;
std::cout << "this should be approx [3.5574 -2.3119]" << std::endl;
// Print the basis and projections: now the new basis
std::cout << std::endl;
std::cout << "Now the basis of projection is: " << std::endl;
for (ImagePCAShapeModelEstimatorType::BasisImagePointerVector::const_iterator
basis_it = basis_check_2.begin(); basis_it != basis_check_2.end(); ++basis_it)
{
std::cout << "[";
InputImageIterator basisImage_it( *basis_it, (*basis_it)->GetBufferedRegion() );
for (basisImage_it.GoToBegin(); !basisImage_it.IsAtEnd(); ++basisImage_it)
{
std::cout << basisImage_it.Get() << " ";
}
std::cout << "]" << std::endl;
}
//Print the projections
std::cout << "The projection of [0 2 2 0] is [" << proj3_2 << "]" << std::endl;
std::cout << "this should be approx [1.4142 -1.4142 0]" << std::endl;
std::cout << "The projection of [0 3 3 0] is [" << proj4_2 << "]" << std::endl;
std::cout << "this should be approx [2.1213 2.1213 3.000]" << std::endl;
//Test for the eigen values for the test case precomputed using Matlab
std::cout << "" << std::endl;
if( proj3[0] < -1.54 && proj3[0] > -1.55 && proj4[1] < -2.31 && proj4[1] > -2.32 &&
proj3_2[1] < -1.414 && proj3_2[1] > -1.415 && proj4_2[2] < 3.01 && proj4_2[2] > 2.99 )
{
std::cerr << "Test Passed" << std::endl;
return EXIT_SUCCESS;
}
else
{
std::cerr << "Test failed" << std::endl;
std::cerr << "The project is out of the range of Matlab precomputed values" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
--Apple-Mail-2-83054363--