[Insight-developers] Re: [Insight-users] One error and one issue with itkImagePCAShapeModelEstimator
Zachary Pincus
zpincus at stanford.edu
Sun, 16 May 2004 04:02:18 -0700
--Apple-Mail-9--976043558
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
> 1) It seems reasonable to allow the ImagePCAShapeModelEstimator
> to accept images of different origin, different StartIndex
> and only enforce the pixel spacing and the image size (in
> pixel along every dimension) to be the same among all the
> input images.
>
> Please send us a patch...
I've included a new implementation of the GenerateInputRequestedRegion
method of itkPCAShapeModelEstimatorClass so images with regions with
different indices, but the same size (and spacing) can work.
Now, the original behavior (which I did not change) is to set all input
requested regions to the largest possible region, and claim to use
that. However, elsewhere in the code (CalculateInnerProducts), the
buffered region of the inputs is what is actually used...
Perhaps it would be better if this were all done with the requested
region? If the user bothered to set them, why blow them away? Or would
that break too much code that relies on not having to set a requested
region?
Anhyow, the diff is pretty messy since this is basically a different
function, so I've just pasted in the new function.
Also attached is a modified test case that exercises this new
functionality, and diffs to the old test.
Zach Pincus
Department of Biochemistry and Program in Biomedical Informatics
Stanford University School of Medicine
New GenerateInputRequestedRegion:
/**
* Requires all of the inputs to have a LargestPosibleRegion and
Spacing the of
* same size.
*/
template<class TInputImage, class TOutputImage>
void
ImagePCAShapeModelEstimator<TInputImage,TOutputImage>
::GenerateInputRequestedRegion()
{
Superclass::GenerateInputRequestedRegion();
// Calling this causes problems because it sets the input requested
regions
// to various output requested regions, which might not overlap,
since the inputs
// and outputs aren't very closely linked in the PCA case.
// If we REALLY want to call this, then some heroics are necessary to
undo
// the damage it does. Specifically, we will have to set all the
inputs'
// requested regions back to *their* LargestPossibleRegion, and not
that of
// some random output.
if ( this->GetInput(0) )
{
// Get the requested region and spacing of the first input
this->GetInput(0)->SetRequestedRegionToLargestPossibleRegion();
typename TInputImage::RegionType masterRequestedRegion =
this->GetInput(0)->GetRequestedRegion();
typename TInputImage::SpacingType masterSpacing =
this->GetInput(0)->GetSpacing();
// Set the requested region of the remaining input to their largest
possible
// regions, and make sure the sizes and spacings are proper.
unsigned int idx;
for (idx = 1; idx < this->GetNumberOfInputs(); ++idx)
{
if ( this->GetInput(idx) )
{
this->GetInput(idx)->SetRequestedRegionToLargestPossibleRegion();
typename TInputImage::RegionType requestedRegion =
this->GetInput(idx)->GetRequestedRegion();
typename TInputImage::SpacingType spacing =
this->GetInput(0)->GetSpacing();
if( masterRequestedRegion.GetSize() !=
requestedRegion.GetSize() ||
masterSpacing != spacing)
{
itkExceptionMacro("Size and spacing of input " << idx <<
" is not the same as that of input 0" );
}
} // if ( this->GetIntput(idx))
} // for idx
} // if( this->GetInput(0) )
}
Diffs for new and old test code
(itkImagePCAShapeModelEstimatorTest.cxx):
< Old
> New
83c83
< InputImageType::IndexType index;
---
> InputImageType::IndexType index, index2;
85c85,87
< InputImageType::RegionType region;
---
> index2.Fill(4);
>
> InputImageType::RegionType region, region2;
89a92,94
> region2.SetSize( inputImageSize );
> region2.SetIndex( index2 );
>
105,106c110,111
< image2->SetLargestPossibleRegion( region );
< image2->SetBufferedRegion( region );
---
> image2->SetLargestPossibleRegion( region2 );
> image2->SetBufferedRegion( region2 );
--Apple-Mail-9--976043558
Content-Transfer-Encoding: 7bit
Content-Type: application/octet-stream;
x-unix-mode=0644;
name="itkImagePCAShapeModelEstimatorTest.cxx"
Content-Disposition: attachment;
filename=itkImagePCAShapeModelEstimatorTest.cxx
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkImagePCAShapeModelEstimatorTest.cxx,v $
Language: C++
Date: $Date: 2004/02/24 09:06:19 $
Version: $Revision: 1.8 $
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 "itkVector.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_math.h"
#include "itkImageRegionIterator.h"
#include "itkLightProcessObject.h"
#include "itkTextOutput.h"
#include "myImagePCAShapeModelEstimator.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::Image<double,NDIMENSION> OutputImageType;
typedef itk::Image<double,NDIMENSION> MeanImageType;
typedef InputImageType::PixelType ImagePixelType;
typedef InputImageType::PixelType InputImagePixelType;
typedef
itk::ImageRegionIterator< InputImageType > InputImageIterator;
typedef
itk::ImageRegionIterator< OutputImageType > OutputImageIterator;
InputImageType::Pointer image1 = InputImageType::New();
InputImageType::Pointer image2 = InputImageType::New();
InputImageType::Pointer image3 = InputImageType::New();
InputImageType::SizeType inputImageSize = {{ IMGWIDTH, IMGHEIGHT }};
InputImageType::IndexType index, index2;
index.Fill(0);
index2.Fill(4);
InputImageType::RegionType region, region2;
region.SetSize( inputImageSize );
region.SetIndex( index );
region2.SetSize( inputImageSize );
region2.SetIndex( index2 );
//--------------------------------------------------------------------------
// Set up Image 1 first
//--------------------------------------------------------------------------
image1->SetLargestPossibleRegion( region );
image1->SetBufferedRegion( region );
image1->Allocate();
// setup the iterators
InputImageIterator image1It( image1, image1->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 2 first
//--------------------------------------------------------------------------
image2->SetLargestPossibleRegion( region2 );
image2->SetBufferedRegion( region2 );
image2->Allocate();
// setup the iterators
InputImageIterator image2It( image2, image2->GetBufferedRegion() );
//--------------------------------------------------------------------------
// Set up Image 3 first
//--------------------------------------------------------------------------
image3->SetLargestPossibleRegion( region );
image3->SetBufferedRegion( region );
image3->Allocate();
// setup the iterators
InputImageIterator image3It( image3, image3->GetBufferedRegion() );
//--------------------------------------------------------------------------
//Manually create and store each vector
//--------------------------------------------------------------------------
//Image no. 1
for( int i = 0; i< 4; i++ )
{
image1It.Set( 1 ); ++image1It;
}
//Image no. 2
image2It.Set( 2 ); ++image2It;
image2It.Set( 0 ); ++image2It;
image2It.Set( 0 ); ++image2It;
image2It.Set( 2 ); ++image2It;
//Image no. 3
image3It.Set( 0 ); ++image3It;
image3It.Set( 3 ); ++image3It;
image3It.Set( 3 ); ++image3It;
image3It.Set( 0 ); ++image3It;
//----------------------------------------------------------------------
// Test code for the Shape model estimator
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//Set the image model estimator
//----------------------------------------------------------------------
typedef itk::ImagePCAShapeModelEstimator<InputImageType, OutputImageType>
ImagePCAShapeModelEstimatorType;
ImagePCAShapeModelEstimatorType::Pointer
applyPCAShapeEstimator = ImagePCAShapeModelEstimatorType::New();
//----------------------------------------------------------------------
//Set the parameters of the clusterer
//----------------------------------------------------------------------
applyPCAShapeEstimator->SetNumberOfTrainingImages( NUMTRAINIMAGES );
applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired( NUMLARGESTPC + 1 );
applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired( NUMLARGESTPC );
applyPCAShapeEstimator->SetInput(0, image1);
applyPCAShapeEstimator->SetInput(1, image2);
applyPCAShapeEstimator->SetInput(2, image3);
applyPCAShapeEstimator->Update();
//Test the printself function to increase coverage
applyPCAShapeEstimator->Print(std::cout);
//Exercise TypeMacro in superclass
typedef ImagePCAShapeModelEstimatorType::Superclass GenericEstimatorType;
std::cout << applyPCAShapeEstimator->GenericEstimatorType::GetNameOfClass() << std::endl;
//Print out the number of training images and the number of principal
//components
std::cout << "The number of training images are: " <<
applyPCAShapeEstimator->GetNumberOfTrainingImages() << std::endl;
std::cout << "The number of principal components desired are: " <<
applyPCAShapeEstimator->GetNumberOfPrincipalComponentsRequired() << std::endl;
//Print the eigen vectors
vnl_vector<double> eigenValues =
applyPCAShapeEstimator->GetEigenValues();
unsigned int numEigVal = eigenValues.size();
std::cout << "Number of returned eign-values: " << numEigVal << std::endl;
std::cout << "The " <<
applyPCAShapeEstimator->GetNumberOfPrincipalComponentsRequired() <<
" largest eigen values are:" << std::endl;
for(unsigned int i= 0; i< vnl_math_min( numEigVal, (unsigned int)NUMLARGESTPC ); i++ )
{
std::cout << eigenValues[ i ] << std::ends;
}
std::cout << "" << std::endl;
std::cout << "" << std::endl;
//Print the MeanImage
OutputImageType::Pointer outImage = applyPCAShapeEstimator->GetOutput( 0 );
OutputImageIterator outImageIt( outImage, outImage->GetBufferedRegion() );
outImageIt.GoToBegin();
std::cout << "The mean image is:" << std::endl;
while(!outImageIt.IsAtEnd() )
{
std::cout << (double)(outImageIt.Get()) << ";" << std::ends;
++outImageIt;
}
std::cout << " " << std::endl;
//Print the largest two eigen vectors
for (unsigned int j=1; j< NUMLARGESTPC + 1; j++ )
{
OutputImageType::Pointer outImage = applyPCAShapeEstimator->GetOutput( j );
OutputImageIterator outImageIt( outImage, outImage->GetBufferedRegion() );
outImageIt.GoToBegin();
std::cout << "" << std::endl;
std::cout << "The eigen vector number: " << j << " is:" << std::endl;
while(!outImageIt.IsAtEnd() )
{
std::cout << (double) (outImageIt.Get()) << ";" << std::ends;
++outImageIt;
}
std::cout << " " << std::endl;
}
//Test for the eigen values for the test case precomputed using Matlab/Splus
std::cout << "" << std::endl;
if( (eigenValues[2] < 6 || eigenValues[2] > 6.1) || (eigenValues[1] >0.1) )
std::cout<< "Test Passed" << std::endl;
else
std::cout<< "Test failed" << std::endl;
return 0;
}
--Apple-Mail-9--976043558
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
charset=US-ASCII;
format=flowed
--Apple-Mail-9--976043558--