[Insight-users] Extract image in a loop, I don't understand why it does not work

[Atwood, Robert C]

 Hi,
I am trying to obtain the statistics for some specified regions of an
image (eventualy each of several images)
 I seem to be getting an exception thrown the second time I try to
extract a region. Is there something that needs to be done to reset the
extractImagefilter? 

In this case I defined the pixel type to unsigned short and the input
image is 1015x512. With this example the size of the region is 20x20 and
the index starts at 0,200 and increases by 5,0 until 20,200 so there
should be lots of room in the image for the region.



Thanks,
Robert

*****Some output ********************

...reading
Region One: ImageRegion (0x7fbfffde70)
  Dimension: 2
  Index: [0, 0]
  Size: [1015, 512]

Region Two: ImageRegion (0x602118)
  Dimension: 2
  Index: [0, 0]
  Size: [1015, 512]

Region Three: ImageRegion (0x7fbfffde70)
  Dimension: 2
  Index: [0, 200]
  Size: [20, 20]

Region Four: ImageRegion (0x602118)
  Dimension: 2
  Index: [0, 0]
  Size: [1015, 512]

Sample mean = [18226.6]
Region Two: ImageRegion (0x602118)
  Dimension: 2
  Index: [0, 0]
  Size: [1015, 512]

Region Three: ImageRegion (0x7fbfffde70)
  Dimension: 2
  Index: [5, 200]
  Size: [20, 20]

Exception caught chopper update !

itk::InvalidRequestedRegionError (0x614500)
Location: "virtual void itk::DataObject::PropagateRequestedRegion()"
File: /sources/local/ITK_cvs/Code/Common/itkDataObject.cxx
Line: 397
Description: Requested region is (at least partially) outside the
largest possible region.


Region Four: ImageRegion (0x602118)
  Dimension: 2
  Index: [0, 0]
  Size: [1015, 512]

Sample mean = [18226.6]




****************** The Code *************************************

#include <stdio.h>

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkScalarImageToListAdaptor.h"
#include "itkScalarToArrayCastImageFilter.h"
#include "itkMeanCalculator.h"
#include "itkExtractImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"

#ifdef T_FLOAT
  typedef float PixelType;
#elif defined T_BYTE
  typedef unsigned char PixelType;
#elif defined T_USHORT
  typedef unsigned short PixelType;
#elif defined T_UINT
  typedef u_int32_t PixelType;
#else
#warning NO TYPE DEFINED, defaulting to FLOAT 
  typedef float PixelType;
#endif


  typedef itk::Image< PixelType, 2 > ImageType;

  // definitions for writing the image
  typedef itk::ImageFileReader<ImageType> ReaderType;
  typedef itk::FixedArray< PixelType, 1 > MeasurementVectorType;
  typedef itk::Image< MeasurementVectorType, 2 > ArrayImageType;
  typedef itk::Statistics::ScalarImageToListAdaptor< ImageType >
SampleType;
  typedef itk::Statistics::MeanCalculator< SampleType >
MeanAlgorithmType;
  typedef itk::ExtractImageFilter<ImageType,ImageType> ChopperType;
/*********************************/
/* beginning of MAIN routine     */
/*********************************/
int main(int argc, char ** argv) {
  int i,j;

  /* ITK objects needed for processing the volume  */
  ReaderType::Pointer reader = ReaderType::New();
  ImageType::Pointer inputImage ;

  ChopperType::Pointer chopper = ChopperType::New();
  SampleType::Pointer sample = SampleType::New();
  MeanAlgorithmType::Pointer meanAlgorithm = MeanAlgorithmType::New();

     ImageType::RegionType inputRegion;
     ImageType::SizeType inputRegionSize;
     ImageType::IndexType inputRegionIndex;

  if (argc != 2){
      printf("Little program to use itk to calculate means of image
area\n");
      printf("Usage: %s input.xxx \n",argv[0]);
      printf("where xxx is an ITK registered file type extension\n");
      return(0);
  }

  /* Read in the  image */
  printf("...reading\n");
  reader->SetFileName(argv[1]);
  inputImage=reader->GetOutput();
  try {
     inputImage->Update();
  }catch( itk::ExceptionObject & exp ) {
     std::cerr << "Exception caught inputImage update !" << std::endl;
     std::cerr << exp << std::endl;
  }
  inputImage->DisconnectPipeline();

  /* select the specified region to get the statistics */
  inputRegion = inputImage->GetLargestPossibleRegion ();
  inputRegionSize = inputRegion.GetSize ();
  inputRegionIndex = inputRegion.GetIndex ();

  /* debugging output */
  std::cout << "Region One: " << inputRegion << std::endl;

  /* loop through regions across the sample */
  chopper->SetInput(inputImage);
  for(i=0;i<20;i+=5){
  std::cout << "**********************************" << std::endl;
  std::cout << "* Region    " << i << "*********************" <<
std::endl;
  std::cout << "**********************************" << std::endl;
     inputRegionIndex[0]=i;
     inputRegionIndex[1]=200;
     inputRegionSize[0]=20;
     inputRegionSize[1]=20;
     inputRegion.SetSize (inputRegionSize);
     inputRegion.SetIndex (inputRegionIndex);

     chopper->SetExtractionRegion (inputRegion);

  /* debugging output */
     std::cout << "Region Two: " <<
inputImage->GetLargestPossibleRegion() << std::endl;
     std::cout << "Region Three: " << inputRegion << std::endl;
     //std::cout << chopper << std::endl;

     try {
        chopper->Update();
     }catch( itk::ExceptionObject & exp ) {
        std::cerr << "Exception caught chopper update !" << std::endl;
        std::cerr << exp << std::endl;
     }
  /* debugging output */
     std::cout << "Region Four: " <<
chopper->GetInput()->GetLargestPossibleRegion() << std::endl;

     /* pipe the image to the statistics package */
     sample->SetImage(chopper->GetOutput());
     //sample->SetImage(inputImage);
     meanAlgorithm->SetInputSample( sample);

     /* get the mean greylevel */
     try {
        meanAlgorithm->Update();
     }catch( itk::ExceptionObject & exp ) {
        std::cerr << "Exception caught subImage update !" << std::endl;
        std::cerr << exp << std::endl;
     }
     std::cout << "Sample mean = " << *(meanAlgorithm->GetOutput()) <<
std::endl;
  }
  printf("All done\n");
  return (0);
}