[ITK] [ITK-users] Picking 10000 samples from a larger mha image using ImageRandomNonRepeatingConstIteratorWithIndex and finding histogram.

[vijayarani]

Hi,
Thanks for your reply.Here is my code and the output

#include "itkImage.h"
#include "itkImageFileReader.h"
#include <itkScalarImageToHistogramGenerator.h>
#include <itkMetaImageIO.h>
#include <itkTIFFImageIO.h>
#include "itkImageFileWriter.h"
#include "itkImageRandomNonRepeatingConstIteratorWithIndex.h"
 
int main(int argc, char*argv[])
{
   const    unsigned int    Dimension = 3;
   typedef unsigned short PixelType;
   typedef itk::Image<PixelType, Dimension>  ImageType;

   typedef itk::ImageFileReader< ImageType  > ImageReaderType;
   typedef itk::ImageFileWriter< ImageType  > ImageWriterType;

   ImageReaderType::Pointer ImageReader   = ImageReaderType::New();
   ImageWriterType::Pointer ImageWriter =  ImageWriterType::New();

   ImageReader->SetFileName( "SampleMha.mha" );
   ImageWriter->SetFileName( "RandomIterator.mha" );

   ImageType::Pointer Image = ImageReader->GetOutput();
   Image->Update();

   ImageType::RegionType region=Image->GetRequestedRegion();
   ImageType::SizeType size=region.GetSize();
   ImageType::IndexType index=region.GetIndex();
   ImageType::IndexType start = index;
   ImageType::SpacingType   spacing   = Image->GetSpacing();
   ImageType::DirectionType direction = Image->GetDirection();

   ImageType::RegionType newregion;
   newregion.SetSize(size);
   newregion.SetIndex(start);
   std::cout << "The size of the image is " << size <<std::endl;
   
   //Define Output Image
   ImageType::Pointer output = ImageType::New();
   output->SetRegions(newregion);
   output->SetSpacing(spacing);
   output->SetDirection(direction);
   output->Allocate();
   output->FillBuffer(0);
   output->Update();
 
   typedef itk::ImageRandomNonRepeatingConstIteratorWithIndex<ImageType >
ConstIteratorType;                                

   ConstIteratorType inputIt(  Image,  region );

   unsigned int number = 10000;
   inputIt.SetNumberOfSamples( number  ); 
   std::cout << "The number of samples is " << number  << std::endl;
   inputIt.ReinitializeSeed();
   inputIt.GoToBegin();

   while(!inputIt.IsAtEnd())  
   {    
     output->SetPixel(inputIt.GetIndex(), inputIt.Get());
     ++inputIt;
   }
   output->Update();
   ImageWriter->SetInput( output );
   ImageWriter->Update();
   calculatePercentile("RandomIterator.mha"); 
   return EXIT_SUCCESS;
}
void calculatePercentile(QString file)
{ 
  qDebug()<<"Calling:  calculatePercentile";
  //Define the pixel type and dimension for reader
  typedef int InternalPixelType;
  const unsigned int    Dimension=3;
  typedef itk::Image<InternalPixelType, Dimension>    InternalImageType;

  //Define the reader
  typedef itk::ImageFileReader<InternalImageType>  ReaderType;
  ReaderType::Pointer reader = ReaderType::New();

  typedef InternalPixelType HistogramMeasurementType;
  typedef itk::Statistics::Histogram< HistogramMeasurementType >
HistogramType;
  typedef itk::Statistics::ImageToListSampleAdaptor< InternalImageType >  
AdaptorType;
  AdaptorType::Pointer adaptor = AdaptorType::New();
  typedef
itk::Statistics::SampleToHistogramFilter<AdaptorType,HistogramType>
FilterType;
  FilterType::Pointer filter = FilterType::New();

  if(QFileInfo(file).suffix()=="mha")
  {
    typedef itk::MetaImageIO    MetaImageIOType;
	MetaImageIOType::Pointer metaIO=MetaImageIOType::New();
	metaIO->SetFileName(file.toLocal8Bit().data());  
	reader->SetFileName(file.toLocal8Bit().data());	
	metaIO->ReadImageInformation();	
	qDebug()<<"No of dimesnsions "<<metaIO->GetNumberOfDimensions()<<"No of
Comps "<<metaIO->GetNumberOfComponents()<<"Comp size
"<<metaIO->GetComponentSize();
    qDebug()<<"Image size in bytes "<<metaIO->GetImageSizeInBytes()<<"Image
size in Comps "<<metaIO->GetImageSizeInComponents()<<"Image size in Pixels
"<<metaIO->GetImageSizeInPixels();
 
	adaptor->SetImage(reader->GetOutput() );
	try
	{
	  reader->Update();
	}
	catch( itk::ExceptionObject & excp )
	{
	  std::cerr << "Problem reading image file" << std::endl;
	  std::cerr << excp << std::endl;  
	}
	const unsigned int numberOfComponents = metaIO->GetNumberOfComponents();
	HistogramType::SizeType size(numberOfComponents);
	size.Fill( metaIO->GetImageSizeInPixels() );	
	filter->SetInput( adaptor );
	filter->SetHistogramSize( size );
	filter->SetMarginalScale( 10 );
	HistogramType::MeasurementVectorType min( numberOfComponents );
	HistogramType::MeasurementVectorType max( numberOfComponents );
	min.Fill( 0 );
	max.Fill( metaIO->GetImageSizeInComponents() );	
	filter->SetHistogramBinMinimum( min );
	filter->SetHistogramBinMaximum( max );
	filter->Update();
  }  
  HistogramType::ConstPointer histogram = filter->GetOutput();
  qDebug()<<"5th percentile is "<< histogram->Quantile( 0, 0.05 ) <<" 25th
percentile is  "<<histogram->Quantile( 0, 0.25 )<<" 50th percentile is 
"<<histogram->Quantile( 0, 0.50 )<<" 75th percentile is 
"<<histogram->Quantile( 0, 0.75 )<<" 95th percentile is
"<<histogram->Quantile( 0, 0.95 );
}

And the output is:
The size of the image is [650, 450, 300]
The number of samples is 10000
Calling:  calculatePercentile
No of dimesnsions  3 No of Comps  1 Comp size  2
Image size in bytes  175500000 Image size in Comps  87750000 Image size in
Pixels  87750000
5th percentile is  0  25th percentile is   0   50th percentile is  0   75th
percentile is  0   95th percentile is 0

I have uploaded the SampleMha.mha file.

Thanks & Regards
Vijaya Rani S



--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/Picking-10000-samples-from-a-larger-mha-image-using-ImageRandomNonRepeatingConstIteratorWithIndex-an-tp7586255p7586259.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users