ITK  4.8.0
Insight Segmentation and Registration Toolkit
WikiExamples/Statistics/ScalarImageKmeansImageFilter.cxx
#include "vtkVersion.h"
#include "vtkImageViewer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkSmartPointer.h"
#include "vtkImageActor.h"
#include "vtkImageMapper3D.h"
#include "vtkInteractorStyleImage.h"
#include "vtkRenderer.h"
static void CreateImage(ImageType::Pointer image);
int main(int, char*[])
{
ImageType::Pointer image = ImageType::New();
CreateImage(image);
KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New();
kmeansFilter->SetInput(image);
kmeansFilter->SetUseNonContiguousLabels(true);
kmeansFilter->AddClassWithInitialMean(8);
kmeansFilter->AddClassWithInitialMean(110);
kmeansFilter->AddClassWithInitialMean(210);
kmeansFilter->Update();
KMeansFilterType::ParametersType estimatedMeans = kmeansFilter->GetFinalMeans();
const unsigned int numberOfClasses = estimatedMeans.Size();
for(unsigned int i = 0 ; i < numberOfClasses ; ++i)
{
std::cout << "cluster[" << i << "] ";
std::cout << " estimated mean : " << estimatedMeans[i] << std::endl;
}
typedef KMeansFilterType::OutputImageType OutputImageType;
OutputImageType,
OutputImageType > RelabelFilterType;
RelabelFilterType::Pointer relabeler = RelabelFilterType::New();
relabeler->SetInput( kmeansFilter->GetOutput() );
RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
rescaleFilter->SetInput(relabeler->GetOutput());
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);
typedef std::vector< unsigned long > SizesType;
const SizesType & sizes = relabeler->GetSizeOfObjectsInPixels();
SizesType::const_iterator sizeItr = sizes.begin();
SizesType::const_iterator sizeEnd = sizes.end();
std::cout << "Number of pixels per class " << std::endl;
unsigned int kclass = 0;
while( sizeItr != sizeEnd )
{
std::cout << "Class " << kclass << " = " << *sizeItr << std::endl;
++kclass;
++sizeItr;
}
// Visualize
ConnectorType::Pointer originalConnector = ConnectorType::New();
originalConnector->SetInput(image);
vtkSmartPointer<vtkImageActor> originalActor =
vtkSmartPointer<vtkImageActor>::New();
#if VTK_MAJOR_VERSION <= 5
originalActor->SetInput(originalConnector->GetOutput());
#else
originalConnector->Update();
originalActor->SetInputData(originalConnector->GetOutput());
#endif
ConnectorType::Pointer outputConnector = ConnectorType::New();
outputConnector->SetInput(rescaleFilter->GetOutput());
vtkSmartPointer<vtkImageActor> outputActor =
vtkSmartPointer<vtkImageActor>::New();
#if VTK_MAJOR_VERSION <= 5
outputActor->SetInput(outputConnector->GetOutput());
#else
outputConnector->Update();
outputActor->SetInputData(outputConnector->GetOutput());
#endif
// There will be one render window
vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();
renderWindow->SetSize(600, 300);
vtkSmartPointer<vtkRenderWindowInteractor> interactor =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
interactor->SetRenderWindow(renderWindow);
// Define viewport ranges
// (xmin, ymin, xmax, ymax)
double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};
// Setup both renderers
vtkSmartPointer<vtkRenderer> leftRenderer =
vtkSmartPointer<vtkRenderer>::New();
renderWindow->AddRenderer(leftRenderer);
leftRenderer->SetViewport(leftViewport);
leftRenderer->SetBackground(.6, .5, .4);
vtkSmartPointer<vtkRenderer> rightRenderer =
vtkSmartPointer<vtkRenderer>::New();
renderWindow->AddRenderer(rightRenderer);
rightRenderer->SetViewport(rightViewport);
rightRenderer->SetBackground(.4, .5, .6);
// Add the sphere to the left and the cube to the right
leftRenderer->AddActor(originalActor);
rightRenderer->AddActor(outputActor);
leftRenderer->ResetCamera();
rightRenderer->ResetCamera();
renderWindow->Render();
vtkSmartPointer<vtkInteractorStyleImage> style =
vtkSmartPointer<vtkInteractorStyleImage>::New();
interactor->SetInteractorStyle(style);
interactor->Start();
return EXIT_SUCCESS;
}
void CreateImage(ImageType::Pointer image)
{
// Create an image with 2 connected components
ImageType::RegionType region;
ImageType::IndexType start;
start[0] = 0;
start[1] = 0;
ImageType::SizeType size;
size[0] = 200;
size[1] = 300;
region.SetSize(size);
region.SetIndex(start);
image->SetRegions(region);
image->Allocate();
itk::ImageRegionIterator<ImageType> imageIterator(image,region);
while(!imageIterator.IsAtEnd())
{
if(imageIterator.GetIndex()[0] > 100 &&
imageIterator.GetIndex()[0] < 150 &&
imageIterator.GetIndex()[1] > 100 &&
imageIterator.GetIndex()[1] < 150)
{
imageIterator.Set(100);
}
else if(imageIterator.GetIndex()[0] > 50 &&
imageIterator.GetIndex()[0] < 70 &&
imageIterator.GetIndex()[1] > 50 &&
imageIterator.GetIndex()[1] < 70)
{
imageIterator.Set(200);
}
else
{
imageIterator.Set(10);
}
++imageIterator;
}
}