[Insight-users] Re: One last call for help - neighborhood operator
dave
common8 at gmail.com
Wed Aug 22 15:19:10 EDT 2007
Skipped content of type multipart/alternative-------------- next part --------------
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkVector.h"
#include "itkListSample.h"
#include "itkKdTree.h"
#include "itkWeightedCentroidKdTreeGenerator.h"
#include "itkKdTreeBasedKmeansEstimator.h"
#include "itkMinimumDecisionRule.h"
#include "itkEuclideanDistance.h"
#include "itkSampleClassifier.h"
#include "itkNormalVariateGenerator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkScalarImageTextureCalculator.h"
#include "itkNeighborhoodIterator.h"
#include "itkRescaleIntensityImageFilter.h"
int main(int argc, char * argv [] )
{
// check input arguments
if( argc < 3 )
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0];
std::cerr << " inputScalarImage outputLabeledScalarImage numberOfClasses(K)" << std::endl;
return EXIT_FAILURE;
}
// parameters
const char * inputImageFileName = argv[1];
const char * outputImageFileName = argv[2];
const unsigned int dimension = 2;
const unsigned int measurementVectorLength = 6;
const unsigned int numberOfClasses = atoi( argv[3] );
// typedefs
typedef char InputPixelType;
typedef unsigned char OutputPixelType;
typedef itk::Image<InputPixelType, dimension > InputImageType;
typedef itk::Image<OutputPixelType, dimension > OutputImageType;
typedef itk::Vector< double,
measurementVectorLength > MeasurementVectorType;
typedef itk::Statistics::ListSample<
MeasurementVectorType > SampleType;
typedef itk::NeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
typedef itk::ImageRegionIterator< OutputImageType > ImageRegionIteratorType;
typedef itk::ImageFileReader< InputImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
typedef itk::Statistics::ScalarImageTextureCalculator<
InputImageType> TexCalcType;
typedef itk::RescaleIntensityImageFilter<
OutputImageType, OutputImageType > RescaleIntensityFilterType;
// Get input image
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( inputImageFileName );
try
{
reader->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// Set up neighborhood iterator and the temp output image
NeighborhoodIteratorType::RadiusType radius;
radius.Fill(1);
NeighborhoodIteratorType it( radius, reader->GetOutput(),
reader->GetOutput()->GetRequestedRegion() );
InputImageType::Pointer tempImage = InputImageType::New();
tempImage->SetRegions( it.GetRegion());
tempImage->Allocate();
// Create texture calculator
TexCalcType::Pointer texCalc = TexCalcType::New();
texCalc->SetInput( tempImage );
texCalc->FastCalculationsOn();
// sample, for storing many measurement vectors
SampleType::Pointer sample = SampleType::New();
sample->SetMeasurementVectorSize( measurementVectorLength );
// Feature Extraction Loop
for (it.GoToBegin() ; !it.IsAtEnd(); ++it)
{
// Feature output information
std::cout << "Calculating features, current center pixel index = " << it.GetIndex() << std::endl;
// Manually update the temp image (this was a last resort because
// the NeighborhoodIterator does not seem to support passing
// the current neighborhood as an image directly)
// This should be fixed for general neighborhood sizes and
// generalized image dimensions!!
NeighborhoodIteratorType::NeighborhoodType neighborhood = it.GetNeighborhood();
InputImageType::IndexType pixelIndex;
pixelIndex[0] = 0;
pixelIndex[1] = 0;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(0) );
pixelIndex[0] = 0;
pixelIndex[1] = 1;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(3) );
pixelIndex[0] = 0;
pixelIndex[1] = 2;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(6) );
pixelIndex[0] = 1;
pixelIndex[1] = 0;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(1) );
pixelIndex[0] = 1;
pixelIndex[1] = 1;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(4) );
pixelIndex[0] = 1;
pixelIndex[1] = 2;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(7) );
pixelIndex[0] = 2;
pixelIndex[1] = 0;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(2) );
pixelIndex[0] = 2;
pixelIndex[1] = 1;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(5) );
pixelIndex[0] = 2;
pixelIndex[1] = 2;
tempImage->SetPixel( pixelIndex, neighborhood.GetElement(8) );
tempImage->Update();
texCalc->Compute();
TexCalcType::FeatureValueVectorPointer means, stds;
means = texCalc->GetFeatureMeans();
// Copy feature means into a measurement vector
// Insert measurement vector into sample (PushBack)
MeasurementVectorType mv;
for( unsigned int i=0; i<measurementVectorLength; ++i)
{
mv[i] = means->GetElement(i);
}
sample->PushBack( mv );
std::cout << " mv = " << mv << std::endl;
}
// Create tree generator
typedef itk::Statistics::WeightedCentroidKdTreeGenerator< SampleType >
TreeGeneratorType;
TreeGeneratorType::Pointer treeGenerator = TreeGeneratorType::New();
treeGenerator->SetSample( sample );
treeGenerator->SetBucketSize( 16 );
treeGenerator->Update();
// Create output tree and kmeans estimator
typedef TreeGeneratorType::KdTreeType TreeType;
typedef itk::Statistics::KdTreeBasedKmeansEstimator<TreeType> EstimatorType;
EstimatorType::Pointer estimator = EstimatorType::New();
// Provide initial feature means
EstimatorType::ParametersType initialMeans(measurementVectorLength*numberOfClasses);
initialMeans.Fill(0.0);
// Set various estimator parameters
estimator->SetParameters( initialMeans );
estimator->SetKdTree( treeGenerator->GetOutput() );
estimator->SetMaximumIteration( 2000 );
estimator->SetCentroidPositionChangesThreshold(0.0);
estimator->StartOptimization();
// Display cluster means
EstimatorType::ParametersType estimatedMeans = estimator->GetParameters();
for ( unsigned int i = 0 ; i < measurementVectorLength ; ++i )
{
std::cout << "cluster[" << i << "] " << std::endl;
std::cout << " estimated mean : " << estimatedMeans[i] << std::endl;
}
// Create membership functions to assign labels from estimator (classifier)
typedef itk::Statistics::EuclideanDistance< MeasurementVectorType >
MembershipFunctionType;
typedef itk::MinimumDecisionRule DecisionRuleType;
DecisionRuleType::Pointer decisionRule = DecisionRuleType::New();
typedef itk::Statistics::SampleClassifier< SampleType > ClassifierType;
ClassifierType::Pointer classifier = ClassifierType::New();
// Set classifier parameters
classifier->SetDecisionRule( (itk::DecisionRuleBase::Pointer) decisionRule);
classifier->SetSample( sample );
classifier->SetNumberOfClasses( numberOfClasses );
// Create class labels (ints from 0 to numberOfClasses)
std::vector< unsigned int > classLabels;
classLabels.resize( numberOfClasses );
for ( unsigned int i = 0; i<numberOfClasses; ++i)
{
classLabels[i] = i;
}
classifier->SetMembershipFunctionClassLabels( classLabels );
// Create a membership function for each classLabel, assign to classifier
std::vector< MembershipFunctionType::Pointer > membershipFunctions;
MembershipFunctionType::OriginType origin( sample->GetMeasurementVectorSize() );
int index = 0;
for ( unsigned int i = 0 ; i < numberOfClasses ; i++ )
{
membershipFunctions.push_back( MembershipFunctionType::New() );
for ( unsigned int j = 0 ; j < sample->GetMeasurementVectorSize(); j++ )
{
origin[j] = estimatedMeans[index++];
}
membershipFunctions[i]->SetOrigin( origin );
classifier->AddMembershipFunction( membershipFunctions[i].GetPointer() );
}
classifier->Update();
// Display measurement vectors and each one's associated classLabel
ClassifierType::OutputType* membershipSample = classifier->GetOutput();
ClassifierType::OutputType::ConstIterator iter = membershipSample->Begin();
while ( iter != membershipSample->End() )
{
std::cout << "measurement vector = " << iter.GetMeasurementVector()
<< "class label = " << iter.GetClassLabel()
<< std::endl;
++iter;
}
// Allocate a 2D image to place the ClassLabels into (same size as input image)
InputImageType::RegionType region = reader->GetOutput()->GetRequestedRegion();
OutputImageType::Pointer classifierOutputImage = OutputImageType::New();
classifierOutputImage->SetRegions( region );
classifierOutputImage->Allocate();
// Set up iterator to extract 1D list of pixels' ClassLabels from membershipSample
ClassifierType::OutputType::ConstIterator iter1 = membershipSample->Begin();
// Set up iterator to place ClassLabels into output image
ImageRegionIteratorType outputIt( classifierOutputImage, region );
outputIt.GoToBegin();
while ( iter1 != membershipSample->End() )
{
outputIt.Set( iter1.GetClassLabel() );
++iter1;
++outputIt;
}
// Rescale intensity of ClassLabels to range of OutputPixelType [char => 0-255]
// This is purely for image display purposes (label #s 1->2 will be very dark)
RescaleIntensityFilterType::Pointer rescaler = RescaleIntensityFilterType::New();
rescaler->SetOutputMinimum( 0 );
rescaler->SetOutputMaximum( 255 );
rescaler->SetInput( classifierOutputImage );
// Write rescaled image to file
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( outputImageFileName );
writer->SetInput( rescaler->GetOutput() );
try
{
writer->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
return 0;
}
More information about the Insight-users
mailing list