ITK/Examples/Statistics/TextureFeatures

From KitwarePublic
< ITK‎ | Examples
Jump to: navigation, search

This example takes a 3D image, and extracts 3 different features for all 13 possible spatial directions. The feature images are not defined along the edges because 3x3x3 window is used to calculate features, and no boundary conditions are employed.

On a Core i7 920 @2.67GHz, it takes about 90 seconds to run this example (with the provided 64x64x12 image) when compiled in release mode using VS2010.

Slice 6 from the original image
Slice 6 from the correlation0 image
Slice 6 from the inertia0 image (with amplified contrast)
Slice 6 from the energy0 image
Slice 6 from the energy12 image

Slices from the original image and 4 images containing derived texture features are displayed on the right.

TextureFeatures.cxx

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <itkDenseFrequencyContainer2.h>
#include "itkHistogramToTextureFeaturesFilter.h"
#include "itkScalarImageToCooccurrenceMatrixFilter.h"
#include "itkVectorContainer.h"
#include "itkAddImageFilter.h"
#include "itkMultiplyImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"

//definitions of used types
typedef itk::Image<float, 3> InternalImageType;
typedef itk::Image<unsigned char, 3> VisualizingImageType;
typedef itk::Neighborhood<float, 3> NeighborhoodType;
typedef itk::Statistics::ScalarImageToCooccurrenceMatrixFilter<InternalImageType>
    Image2CoOccuranceType;
typedef Image2CoOccuranceType::HistogramType HistogramType;
typedef itk::Statistics::HistogramToTextureFeaturesFilter<HistogramType> Hist2FeaturesType;
typedef InternalImageType::OffsetType OffsetType;
typedef itk::AddImageFilter <InternalImageType> AddImageFilterType;
typedef itk::MultiplyImageFilter<InternalImageType> MultiplyImageFilterType;

//calculate features for one offset
void calcTextureFeatureImage (OffsetType offset,
    InternalImageType::Pointer inputImage, InternalImageType::Pointer outInertia,
    InternalImageType::Pointer outCorrelation, InternalImageType::Pointer outEnergy)
{
    //allocate output images
    outInertia->CopyInformation(inputImage);
    outInertia->SetRegions(inputImage->GetLargestPossibleRegion());
    outInertia->Allocate();
    outInertia->FillBuffer(0);
    outCorrelation->CopyInformation(inputImage);
    outCorrelation->SetRegions(inputImage->GetLargestPossibleRegion());
    outCorrelation->Allocate();
    outCorrelation->FillBuffer(0);
    outEnergy->CopyInformation(inputImage);
    outEnergy->SetRegions(inputImage->GetLargestPossibleRegion());
    outEnergy->Allocate();
    outEnergy->FillBuffer(0);

    Image2CoOccuranceType::Pointer glcmGenerator=Image2CoOccuranceType::New();
    glcmGenerator->SetOffset(offset);
    glcmGenerator->SetNumberOfBinsPerAxis(16); //reasonable number of bins
    glcmGenerator->SetPixelValueMinMax(0, 255); //for input UCHAR pixel type
    Hist2FeaturesType::Pointer featureCalc=Hist2FeaturesType::New();

    typedef itk::RegionOfInterestImageFilter<InternalImageType,InternalImageType> roiType;
    roiType::Pointer roi=roiType::New();
    roi->SetInput(inputImage);

    InternalImageType::RegionType window;
    InternalImageType::RegionType::SizeType size;
    size.Fill(3); //window size=3x3x3
    window.SetSize(size);
    InternalImageType::IndexType pi; //pixel index
    
    //slide window over the entire image
    for (unsigned x=1; x<inputImage->GetLargestPossibleRegion().GetSize(0)-1; x++)
    {
        pi.SetElement(0,x);
        window.SetIndex(0,x-1);
        for (unsigned y=1; y<inputImage->GetLargestPossibleRegion().GetSize(1)-1; y++)
        {
            pi.SetElement(1,y);
            window.SetIndex(1,y-1);
            for (unsigned z=1; z<inputImage->GetLargestPossibleRegion().GetSize(2)-1; z++)
            {
                pi.SetElement(2,z);
                window.SetIndex(2,z-1);
                roi->SetRegionOfInterest(window);
                roi->Update();
                glcmGenerator->SetInput(roi->GetOutput());
                glcmGenerator->Update();
                featureCalc->SetInput( glcmGenerator->GetOutput() );
                featureCalc->Update();

                outInertia->SetPixel(pi, featureCalc->GetFeature(Hist2FeaturesType::Inertia));
                outCorrelation->SetPixel(pi, featureCalc->GetFeature(Hist2FeaturesType::Correlation));
                outEnergy->SetPixel(pi, featureCalc->GetFeature(Hist2FeaturesType::Energy));
            }
        }
        std::cout<<'.';
    }
}

int main(int argc, char*argv[])
{
  if(argc < 2)
    {
    std::cerr << "Usage: " << argv[0] << " Required image.mha" << std::endl;
    return EXIT_FAILURE;
    }
  
  std::string fileName = argv[1];
  
  typedef itk::ImageFileReader<InternalImageType> ReaderType;
  ReaderType::Pointer reader=ReaderType::New();
  reader->SetFileName(fileName);
  reader->Update();
  InternalImageType::Pointer image=reader->GetOutput();

  NeighborhoodType neighborhood;
  neighborhood.SetRadius(1);
  unsigned int centerIndex = neighborhood.GetCenterNeighborhoodIndex();
  OffsetType offset;

  typedef itk::ImageFileWriter<InternalImageType> WriterType;
  WriterType::Pointer writer=WriterType::New();

  for ( unsigned int d = 0; d < centerIndex; d++ )
  {
      offset = neighborhood.GetOffset(d);
      InternalImageType::Pointer inertia=InternalImageType::New();
      InternalImageType::Pointer correlation=InternalImageType::New();
      InternalImageType::Pointer energy=InternalImageType::New();
      calcTextureFeatureImage(offset, image, inertia, correlation, energy);
      
      writer->SetInput(inertia);
//	snprintf(buf, 100, "Inertia%u.mha", d); // Warning: call to int __builtin___snprintf_chk will always overflow destination buffer
      std::stringstream ssInertia;
      ssInertia << "Inertia" << d << ".mha";
      writer->SetFileName(ssInertia.str());
      writer->Update();
      writer->SetInput(correlation);
      std::stringstream ssCorrelation;
      ssCorrelation << "Correlation" << d << ".mha";
      writer->SetFileName(ssCorrelation.str());
      writer->Update();
      writer->SetInput(energy);
      std::stringstream ssEnergy;
      ssEnergy << "Energy" << d << ".mha";
      writer->SetFileName(ssEnergy.str());
      writer->Update();
      std::cout<<'\n';
  }

  return EXIT_SUCCESS;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.9.5)

project(TextureFeatures)

find_package(ITK REQUIRED)
include(${ITK_USE_FILE})
if (ITKVtkGlue_LOADED)
  find_package(VTK REQUIRED)
  include(${VTK_USE_FILE})
endif()

add_executable(TextureFeatures MACOSX_BUNDLE TextureFeatures.cxx)

if( "${ITK_VERSION_MAJOR}" LESS 4 )
  target_link_libraries(TextureFeatures ITKReview ${ITK_LIBRARIES})
else( "${ITK_VERSION_MAJOR}" LESS 4 )
  target_link_libraries(TextureFeatures ${ITK_LIBRARIES})
endif( "${ITK_VERSION_MAJOR}" LESS 4 )

Download and Build TextureFeatures

Click here to download TextureFeatures and its CMakeLists.txt file. Once the tarball TextureFeatures.tar has been downloaded and extracted,

cd TextureFeatures/build
  • If ITK is installed:
cmake ..
  • If ITK is not installed but compiled on your system, you will need to specify the path to your ITK build:
cmake -DITK_DIR:PATH=/home/me/itk_build ..

Build the project:

make

and run it:

./TextureFeatures

WINDOWS USERS PLEASE NOTE: Be sure to add the ITK bin directory to your path. This will resolve the ITK dll's at run time.

Building All of the Examples

Many of the examples in the ITK Wiki Examples Collection require VTK. You can build all of the the examples by following these instructions. If you are a new VTK user, you may want to try the Superbuild which will build a proper ITK and VTK.

ItkVtkGlue

ITK >= 4

For examples that use QuickView (which depends on VTK), you must have built ITK with Module_ITKVtkGlue=ON.

ITK < 4

Some of the ITK Examples require VTK to display the images. If you download the entire ITK Wiki Examples Collection, the ItkVtkGlue directory will be included and configured. If you wish to just build a few examples, then you will need to download ItkVtkGlue and build it. When you run cmake it will ask you to specify the location of the ItkVtkGlue binary directory.