ITK/Examples/Statistics/TextureFeatures

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.

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

TextureFeatures.cxx

<source lang="cpp">

1. include "itkImage.h"
2. include "itkImageFileReader.h"
3. include "itkImageFileWriter.h"
4. include <itkDenseFrequencyContainer2.h>
5. include "itkHistogramToTextureFeaturesFilter.h"
6. include "itkScalarImageToCooccurrenceMatrixFilter.h"
7. include "itkVectorContainer.h"
8. include "itkAddImageFilter.h"
9. include "itkMultiplyImageFilter.h"
10. 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;

} </source>

CMakeLists.txt

<syntaxhighlight lang="cmake"> 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 )

</syntaxhighlight>

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.