[Insight-developers] Submission: PCA decomposition calculator
Zachary Pincus
zpincus at stanford.edu
Thu, 29 Apr 2004 18:45:59 -0700
--Apple-Mail-9--244339601
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
I have written a calculator class that projects images down to a space
spanned by a specified orthonormal basis set and returns that
projection. (Also test code too!)
This task is generically useful, but it is especially used in image
processing for computing the "PCA decomposition coefficients" of
images, given a set of some number of principal components. (As a
reminder, the principal components are an orthonormal basis set
spanning some subspace of the original image space, the "decomposition
coefficients" vector is the projection of the image down to this
subspace.)
People tend to use PCA decompositions for image classification tasks,
because a decomposition captures how an image varies along "important"
dimensions, and disregards noise. (cf. the "eigenfaces" literature for
face recognition.)
Hopefully this code will prove useful and merits inclusion in ITK.
Zach Pincus
Department of Biochemistry and Program in Biomedical Informatics
Stanford University School of Medicine
--Apple-Mail-9--244339601
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
x-unix-mode=0644;
name="itkImagePCADecompositionCalculator.h"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculator.h
#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. This is not
* checked at runtime.
*
* \ingroup Operators
*/
template <class TInputImage,
class TBasisImage = Image<double, TInputImage::ImageDimension> >//::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-9--244339601
Content-Transfer-Encoding: 7bit
Content-Type: application/octet-stream;
x-unix-mode=0644;
name="itkImagePCADecompositionCalculator.txx"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculator.txx
#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 (this->m_BasisImages != _arg)
{
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) CalculateBasisMatrix();
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;
}
/*
* 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!");
}
m_ImageAsVector = BasisVectorType(m_NumPixels);
ImageRegionConstIterator<InputImageType> image_it(m_Image,
m_Image->GetRequestedRegion());
int i = 0;
for (image_it.GoToBegin(); !image_it.IsAtEnd(); ++image_it)
{
m_ImageAsVector(i++) = 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-9--244339601
Content-Transfer-Encoding: 7bit
Content-Type: application/octet-stream;
x-unix-mode=0644;
name="itkImagePCADecompositionCalculatorTest.cxx"
Content-Disposition: attachment;
filename=itkImagePCADecompositionCalculatorTest.cxx
#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"
//Data definitions
#define IMGWIDTH 2
#define IMGHEIGHT 2
#define NDIMENSION 2
#define NUMTRAINIMAGES 3
#define NUMLARGESTPC 2
// 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());
//------------------------------------------------------
//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::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() );
//--------------------------------------------------------------------------
//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
//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
image3It.Set( 0 ); ++image3It;
image3It.Set( 3 ); ++image3It;
image3It.Set( 3 ); ++image3It;
image3It.Set( 0 ); ++image3It;
//----------------------------------------------------------------------
// 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();
//Test the printself function to increase coverage
decomposer->Print(std::cout);
// Retrieve the basis images
std::cout << "The basis of projection is: " << std::endl;
ImagePCAShapeModelEstimatorType::BasisImagePointerVector basis2;
basis2 = decomposer->GetBasisImages();
for (ImagePCAShapeModelEstimatorType::BasisImagePointerVector::const_iterator
basis_it = basis2.begin(); basis_it != basis2.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;
//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 )
{
std::cout<< "Test Passed" << std::endl;
return 0;
}
else
{
std::cout<< "Test failed" << std::endl;
return -1;
}
return 0;
}
--Apple-Mail-9--244339601--