[Insight-users] One error and one issue with itkImagePCAShapeModelEstimator

[Zachary Pincus]

--Apple-Mail-3-873457500
Content-Transfer-Encoding: 7bit
Content-Type: text/plain;
	charset=US-ASCII;
	delsp=yes;
	format=flowed

Hello,

I've run into a problematic situation: I want to estimate a PCA model  
of several regions (same size) from the same image. I just run an  
ExtractImageFilter a few times to get a bunch of sub-images, and feed  
them to the model estimator. This should be no problem, but:

(1) itk::ImagePCAShapeModelEstimator will not accept images with  
different LargestPossibleRegions, regardless of whether they are the  
same size. This seems unnecessarily restrictive to me. If others agree,  
I'll happily submit a patch to remedy this behavior.

The more problematic thing is that I had a really hard time tracking  
down this bug because:
(2) Somehow the exception reporting from this class seems to be broken.  
I'm not sure if it's just me, but even though I've got a try/catch  
block, I can't catch the exception raised by an  
ImagePCAShapeModelEstimator object when you try to feed it wrong-sized  
images. This exception is raised on line 146 of  
itkImagePCAShapeModelEstimator.txx -- I know this because when I put a  
bunch of printfs (well, std::cout <<, but you get my drift) in the  
template code, I get output right up until where the exception is to be  
called. Somehow, the exception doesn't get propagated back properly.

I've included some stripped-down test code (derived from the test code  
in the ITK tree) that on my machine reproduces this problem. The code  
tries to feed an ImagePCAShapeModelEstimator images with different  
regions (though the sizes are the same). All of this is done within a  
try/catch block; however the exception raised on line 146 of  
itkImagePCAShapeModelEstimator.txx never gets caught. (My setup: OS X  
10.3.3, GCC 3.3 using c++ (not g++) compiler, with a CVS version of ITK  
updated last night, built with shared libraries.)

Strange,

Zach Pincus


Attached is the test code. Here follows a backtrace of what I get when  
I run it in GDB. There is some strangeness afoot even in the backtrace  
-- note the "itkOctreeNode" source files erroneously assigned to  
symbols from the standard c++ library.

#0  0x90042aac in kill ()
#1  0x9009ec5c in abort ()
#2  0x00018034 in __cxxabiv1::__terminate(void (*)()) () at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkSmartPointer.h:66
#3  0x00018078 in std::terminate() () at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkSmartPointer.h:66
#4  0x028079fc in __cxxabiv1::__unexpected(void (*)()) () at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkOctreeNode.h:117
#5  0x027f607c in __cxa_call_unexpected () at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkOctreeNode.h:117
#6  0x027d3e94 in itk::DataObject::PropagateRequestedRegion()  
(this=0x2501110) at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkDataObject.h:90
#7  0x027d3388 in itk::DataObject::Update() (this=0x25009e0) at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkDataObject.cxx:343
#8  0x0000343c in main () at  
/Volumes/Atlotl/Developer/ITK/ITK-CVS-src/Insight/Code/Common/ 
itkSmartPointer.h:71


--Apple-Mail-3-873457500
Content-Transfer-Encoding: 7bit
Content-Type: application/text;
	x-mac-type=54455854;
	x-unix-mode=0644;
	name="itkImagePCAShapeModelEstimatorTest.cxx"
Content-Disposition: attachment;
	filename=itkImagePCAShapeModelEstimatorTest.cxx

// Insight classes
#include "itkImage.h"
#include "itkImageRegionIterator.h"

#include "itkImagePCAShapeModelEstimator.h"

//Data definitions 
#define   IMGWIDTH            2
#define   IMGHEIGHT           2
#define   NDIMENSION          2
#define   NUMTRAINIMAGES      3
#define   NUMLARGESTPC        2


int main(int, char* [] )
{

  //------------------------------------------------------
  //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(100);
  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( region );
  image2->SetBufferedRegion( region );
  image2->Allocate();

  // setup the iterators
  InputImageIterator image2It( image2, image2->GetBufferedRegion() );

  //--------------------------------------------------------------------------
  // Set up Image 3 first
  //--------------------------------------------------------------------------

  image3->SetLargestPossibleRegion( region2 );
  image3->SetBufferedRegion( region2 );
  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
  //----------------------------------------------------------------------
  try {
	  
  applyPCAShapeEstimator->SetNumberOfTrainingImages( NUMTRAINIMAGES );
  applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired( NUMLARGESTPC );
  
  applyPCAShapeEstimator->SetInput(0, image1);
  applyPCAShapeEstimator->SetInput(1, image2);
  applyPCAShapeEstimator->SetInput(2, image3);
 
  applyPCAShapeEstimator->Update();

  } catch( itk::ExceptionObject & err ) { 
	  std::cout << "ExceptionObject caught !" << std::endl; 
	  std::cout << err << std::endl; 
	  return -1;
  }

  return 0;
}

--Apple-Mail-3-873457500--