[ITK] Grey Level Size Zone Matrix
Cyril Jaudet
drcjaudet at gmail.com
Sat Apr 22 04:42:49 EDT 2017
So, here the code.
You need an image and a label. The calculation of the texture will be into
the label. You normally have to preprocess the image to a specific number
of greylevel, usualy 256 for CT, 64 for PET...
They can be adjusted with maxValue and minValue. MaxVolume is the maximal
volume in number of voxel that you expect a zone will have, if you don't
now just use the volume of the lesion but it can be process consuming.
Also i crop the image and the label before starting the analyse to spare
some ressource. It's a small part f a big program so i hope i didn't forget
any part.
Tell me if you have problem with it or see an error in the code.
Best regards,
Cyril Jaudet, PhD
Medical Physicist
UZBrussel, Radiotherapy
You will have to initialise all the filter but that itk basic stuff like:
°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
#include "itkBinaryThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkLabelGeometryImageFilter.h"
//maybe i have forgetted some just add then in this case.
typedef itk::BinaryThresholdImageFilter <ImageType, LabelType>
BinaryThresholdImageFilterType2_GSZM;
typedef itk::ConnectedComponentImageFilter <LabelType, LabelType>
ConnectedComponentImageFilterType_GSZM;
typedef itk::LabelGeometryImageFilter<LabelType,ImageType >
LabelGeometryImageFilterType_GSZM;
//biblio general
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <dirent.h>
#include <unistd.h>
#include <errno.h>
#include <string.h>
#include <cstdlib>
#include <list>
#include <fstream>
#include <sys/types.h>
#include <sys/stat.h>
#include <algorithm>
#include <math.h>
#include <vector>
const unsigned int ImageDimension = 3;
typedef float PixelType;
typedef itk::Image<PixelType, ImageDimension> ImageType;
typedef itk::Image<unsigned int, ImageDimension> LabelType;
//initialise les vecteurs:
int num_of_col = Max_Value+1;
int num_of_row = Max_Volume+1;
int init_value = 0;
// std::cout << " num_of_col: " << num_of_col << std::endl;
// std::cout << " num_of_row: " << num_of_row << std::endl;
// std::cout << " Rmax_init: " << Rmax_init << std::endl;
vector<int> volume_all ;
//now we have an empty 1D-matrix of size (0). Resizing it with one
single command:
volume_all.resize( num_of_col, 0 );
vector< vector<int> > ZM;
//now we have an empty 2D-matrix of size (0,0). Resizing it with one
single command:
ZM.resize( num_of_col , vector<int>( num_of_row , init_value ) );
//fill the matrix for each distance with the number of grey value
vector< vector<int> > matrix_GLDZM;
//now we have an empty 2D-matrix of size (0,0). Resizing it with one
single command:
matrix_GLDZM.resize( Max_Value+1 , vector<int>(Rmax_init+1 ,
init_value) );
//----------------------------------------------GSZM:
ImageType::Pointer image_GSZM = cropFilter->GetOutput();
image_GSZM->Update();
LabelType::Pointer label_GSZM = LabelcropFilter->GetOutput();
label_GSZM->Update();
int MaxValue=256;
int MinValue=0;
//////////////////////////////////////////:::Connected composants
analysis///////////////////////////////////:
// intiate filter:
BinaryThresholdImageFilterType2_GSZM::Pointer
thresholdFilterToBinary_GSZM= BinaryThresholdImageFilterType2_GSZM::New();
ConnectedComponentImageFilterType_GSZM::Pointer connected_GSZM =
ConnectedComponentImageFilterType_GSZM::New ();
LabelGeometryImageFilterType_GSZM::Pointer labelGeometryImageFilter_GSZM =
LabelGeometryImageFilterType_GSZM::New();
//Maximum volume of the mask:
int D=int(maxValue-minValue);
//std::cout << " Print D : " << D << std::endl;
// int volume_all[int(D)];
float R;
//now we have an empty 1D-matrix of size (0). Refill it with one single
command:
std::fill(volume_all.begin(), volume_all.end(), 0);
/*for(int k=minValue; k<(maxValue+1); k++){ //all resampled Value are
between 0 and D
volume_all[k]=0;
std::cout << " Initiale levels : " << k <<" volume : "<<volume_all[k]
<< std::endl;
}*/
itk::ImageRegionIterator<ImageType>ItI_GSZM(image_GSZM,image_GSZM->GetLargestPossibleRegion());
itk::ImageRegionIterator<LabelType>ItM_GSZM(label_GSZM,label_GSZM->GetLargestPossibleRegion());
int level=0;
for ( ItI_GSZM.GoToBegin(), ItM_GSZM.GoToBegin();!ItI_GSZM.IsAtEnd();
++ItI_GSZM, ++ItM_GSZM )
{
if ( int(ItM_GSZM.Get()) ==A_label )
{
level=ItI_GSZM.Get();
volume_all[level]=volume_all[level]+1;
}
}
int mask_volume=0;
int level_mask_volume=0;
for(int k=minValue; k<(maxValue+1); k++){
if (mask_volume<volume_all[k])
{
mask_volume=volume_all[k];
level_mask_volume=k;
}
// std::cout << " levels : " << k <<" volume : "<<volume_all[k] <<
std::endl;
}
//std::cout << " volume mask : " << mask_volume <<" voxels"<< std::endl;
//std::cout << " level volume mask : " << level_mask_volume << std::endl;
//now we have an empty 2D-matrix of size (0,0). Refill it with one single
command:
ZM.assign( num_of_col, vector<int>(num_of_row) ) ;
//std::cout << "MaskImage initialaise 1"<< std::endl;
int volume=0;
int Max_volume=0;
int Min_volume=mask_volume;
int N_object_connect=0;
for(int level_Grey=minValue; level_Grey<(maxValue+1); level_Grey++ )
{
//choose one levels of Grey (D levels)
thresholdFilterToBinary_GSZM->SetInput(image_GSZM);
thresholdFilterToBinary_GSZM->SetLowerThreshold(int(level_Grey));
thresholdFilterToBinary_GSZM->SetUpperThreshold(int(level_Grey));
thresholdFilterToBinary_GSZM->SetInsideValue(int(level_Grey));
thresholdFilterToBinary_GSZM->SetOutsideValue(0);
thresholdFilterToBinary_GSZM->Update();
//look at the connectivity between voxels each zone is save in a different
label
connected_GSZM->SetInput(thresholdFilterToBinary_GSZM->GetOutput());
connected_GSZM->SetFullyConnected(1);
connected_GSZM->Update();
//std::cout << "Number of objects for D= "<<level_Grey <<": "<<
connected_GSZM->GetObjectCount() << std::endl;
N_object_connect=0;
N_object_connect=connected_GSZM->GetObjectCount();
//measure the volume of each labels
labelGeometryImageFilter_GSZM->SetInput(connected_GSZM->GetOutput() );
labelGeometryImageFilter_GSZM->SetIntensityInput(image_GSZM);
labelGeometryImageFilter_GSZM->Update();
//store in the Thibault matrix
for(int N=1; N<int(N_object_connect+1);N++)
{
volume=labelGeometryImageFilter_GSZM->GetVolume(N);
ZM[level_Grey][volume]++;
// std::cout << "Grey value of: "<<level_Grey+1<< " Matrix filling
volume: " << volume <<" voxels"<< std::endl;
if ( Max_volume<volume ) {Max_volume=volume;}
if ( Min_volume>volume ) {Min_volume=volume;}
}
//inputImage->DisconnectPipeline();
}
//std::cout << "GSZM create"<< std::endl;
std::cout << "max volume :"<< Max_volume<< std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////:
/////////////////////convert size zone matrix to an histogram//////
typedef itk::Statistics::SparseFrequencyContainer2
FrequencyContainerType_GSZM;
typedef FrequencyContainerType_GSZM::AbsoluteFrequencyType
FrequencyType_GSZM;
typedef itk::Statistics::Histogram< float, FrequencyContainerType_GSZM>
HistogramType_GSZM;
unsigned int numberOfComponents=2;
HistogramType_GSZM::Pointer histogram = HistogramType_GSZM::New();
HistogramType_GSZM::SizeType size(numberOfComponents);
size[0]=D;
size[1]=(Max_volume-Min_volume);
HistogramType_GSZM::MeasurementVectorType
lowerBound(numberOfComponents);
lowerBound[0]=0;
lowerBound[1]=0;
HistogramType_GSZM::MeasurementVectorType
upperBound(numberOfComponents);
upperBound[0]=(maxValue-minValue);
upperBound[1]=Max_volume-1;
histogram->SetMeasurementVectorSize(numberOfComponents);
histogram->Initialize(size, lowerBound, upperBound );
histogram->SetToZero();
// std::cout << "Histogramme created: middle "<< std::endl;
HistogramType_GSZM::IndexType index(numberOfComponents);
for(int level_Grey=minValue; level_Grey<(maxValue+1); level_Grey++ )
{
for(int N=0; N<(Max_volume);N++)
{
index[0]=(level_Grey-minValue); //the calulation add +1 to the
two index
index[1]=N;
histogram->SetFrequencyOfIndex(index,
static_cast<FrequencyType_GSZM>(ZM[level_Grey][N+1]));
/* if (N==(1) )
{
std::cout << "index: "<<index[0]<<" "<<index[1]<<" value: "<<
Z[level_Grey][N]<< std::endl;
}*/
}
}
// std::cout << "Histogramme created: ok "<< std::endl;
//-------------------libére les vecteurs:
//ZM.clear();
//////////////////////////////////calcul of matrix size
features////////////////
/////////////:based on the run lenght matrix////////////////:
typedef
itk::Statistics::HistogramToRunLengthFeaturesFilter<HistogramType_GSZM>
HistoSizeFilterType;
HistoSizeFilterType::Pointer hzf= HistoSizeFilterType::New();
hzf->SetInput(histogram);
//hzf->FastCalculationsOn();//just one direction
hzf->Update();
std::cout<<" GSZM: ok"<<std::endl;
//And then to acess the values
//GSZM:
fichier<<hzf->GetShortRunEmphasis()<<"\t";
fichier<<hzf->GetLongRunEmphasis()<<"\t";
fichier<<hzf->GetGreyLevelNonuniformity()<<"\t";
fichier<<hzf->GetRunLengthNonuniformity()<<"\t";
fichier<<hzf->GetLowGreyLevelRunEmphasis()<<"\t";
fichier<<hzf->GetHighGreyLevelRunEmphasis()<<"\t";
fichier<<hzf->GetShortRunLowGreyLevelEmphasis()<<"\t";
fichier<<hzf->GetShortRunHighGreyLevelEmphasis()<<"\t";
fichier<<hzf->GetLongRunLowGreyLevelEmphasis()<<"\t";
fichier<<hzf->GetLongRunHighGreyLevelEmphasis()<<"\t";
°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
2017-04-21 19:31 GMT+02:00 Pati, Sarthak <Sarthak.Pati at uphs.upenn.edu>:
> Hi Cyril,
>
>
>
> Sure, that would be fine, as long as there is nothing missing.
>
>
>
> Thanks so much,
>
> Sarthak
>
>
>
> *From:* Cyril Jaudet [mailto:drcjaudet at gmail.com]
> *Sent:* 21/Apr/2017; Friday 13:31
> *To:* Pati, Sarthak <Sarthak.Pati at uphs.upenn.edu>
> *Subject:* Re: Grey Level Size Zone Matrix
>
>
>
> Hello Pati,
>
> I just have the code that work fine. I'll send it to you tomorrow. It is
> ok if i just give the corresponding part because now it is a part of a
> bigger code.
>
> See you,
>
> Cyril Jaudet, PhD
>
> Medical physicist
>
> UZBrussel, Belgium
>
>
>
> 2017-04-21 18:10 GMT+02:00 Pati, Sarthak <Sarthak.Pati at uphs.upenn.edu>:
>
> Hi,
>
>
>
> I saw your post on the ITK forums (http://public.kitware.com/
> pipermail/community/2015-January/008102.html) and I wanted to know if you
> had that code available somewhere? I was looking for an implementation for
> a project and would gladly refer to your work for this.
>
>
>
> Thanks,
>
> Sarthak Pati
>
> Lead Software Developer
>
> Center for Biomedical Image Computing & Analytics
>
> University of Pennsylvania
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20170422/487e0c05/attachment-0001.html>
More information about the Community
mailing list