<div dir="ltr"><div><div>So, here the code.<br></div>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...<br></div><div>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.<br></div><div>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.<br></div><div><br><div>Tell me if you have problem with it or see an error in the code.<br></div><div>Best regards, <br></div><div>Cyril Jaudet, PhD<br></div><div>Medical Physicist<br></div>UZBrussel, Radiotherapy<br><br><br><br>You will have to initialise all the filter but that itk basic stuff like:<br><br>°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°<br></div><div>#include "itkBinaryThresholdImageFilter.h"<br>#include "itkConnectedComponentImageFilter.h"<br>#include "itkLabelGeometryImageFilter.h"<br><br>//maybe i have forgetted some just add then in this case.<br><br><br>typedef itk::BinaryThresholdImageFilter <ImageType, LabelType> BinaryThresholdImageFilterType2_GSZM;<br><br>typedef itk::ConnectedComponentImageFilter <LabelType, LabelType> ConnectedComponentImageFilterType_GSZM;<br><br>typedef itk::LabelGeometryImageFilter<LabelType,ImageType > LabelGeometryImageFilterType_GSZM;<br><br>//biblio general<br>#include <stdio.h><br>#include <stdlib.h><br>#include <string.h><br>#include <sys/types.h><br>#include <dirent.h><br>#include <unistd.h><br>#include <errno.h><br>#include <string.h><br>#include <cstdlib><br>#include <list><br>#include <fstream><br>#include <sys/types.h><br>#include <sys/stat.h><br>#include <algorithm><br>#include <math.h><br>#include <vector><br><br>const unsigned int ImageDimension = 3;<br> typedef float PixelType;<br><br> typedef itk::Image<PixelType, ImageDimension> ImageType;<br><br> typedef itk::Image<unsigned int, ImageDimension> LabelType;<br><br><br><br> //initialise les vecteurs:<br><br> int num_of_col = Max_Value+1;<br> int num_of_row = Max_Volume+1;<br> int init_value = 0;<br><br> // std::cout << " num_of_col: " << num_of_col << std::endl;<br> // std::cout << " num_of_row: " << num_of_row << std::endl;<br> // std::cout << " Rmax_init: " << Rmax_init << std::endl;<br><br><br><br> vector<int> volume_all ;<br> //now we have an empty 1D-matrix of size (0). Resizing it with one single command:<br> volume_all.resize( num_of_col, 0 );<br> vector< vector<int> > ZM;<br> //now we have an empty 2D-matrix of size (0,0). Resizing it with one single command:<br> ZM.resize( num_of_col , vector<int>( num_of_row , init_value ) );<br><br> //fill the matrix for each distance with the number of grey value<br> vector< vector<int> > matrix_GLDZM;<br> //now we have an empty 2D-matrix of size (0,0). Resizing it with one single command:<br> matrix_GLDZM.resize( Max_Value+1 , vector<int>(Rmax_init+1 , init_value) );<br><br><br><br>//----------------------------------------------GSZM:<br><br><br> ImageType::Pointer image_GSZM = cropFilter->GetOutput();<br> image_GSZM->Update();<br><br> LabelType::Pointer label_GSZM = LabelcropFilter->GetOutput();<br> label_GSZM->Update();<br><br></div><div>int MaxValue=256;<br></div><div>int MinValue=0;<br></div><div><br>//////////////////////////////////////////:::Connected composants analysis///////////////////////////////////:<br><br>// intiate filter:<br><br> BinaryThresholdImageFilterType2_GSZM::Pointer thresholdFilterToBinary_GSZM= BinaryThresholdImageFilterType2_GSZM::New();<br> ConnectedComponentImageFilterType_GSZM::Pointer connected_GSZM = ConnectedComponentImageFilterType_GSZM::New ();<br> LabelGeometryImageFilterType_GSZM::Pointer labelGeometryImageFilter_GSZM = LabelGeometryImageFilterType_GSZM::New();<br><br> //Maximum volume of the mask:<br><br> int D=int(maxValue-minValue);<br>//std::cout << " Print D : " << D << std::endl;<br><br>// int volume_all[int(D)];<br> float R;<br><br> //now we have an empty 1D-matrix of size (0). Refill it with one single command:<br><br> std::fill(volume_all.begin(), volume_all.end(), 0);<br><br> /*for(int k=minValue; k<(maxValue+1); k++){ //all resampled Value are between 0 and D<br> volume_all[k]=0;<br> std::cout << " Initiale levels : " << k <<" volume : "<<volume_all[k] << std::endl;<br> }*/<br><br><br> itk::ImageRegionIterator<ImageType>ItI_GSZM(image_GSZM,image_GSZM->GetLargestPossibleRegion());<br> itk::ImageRegionIterator<LabelType>ItM_GSZM(label_GSZM,label_GSZM->GetLargestPossibleRegion());<br> int level=0;<br><br> for ( ItI_GSZM.GoToBegin(), ItM_GSZM.GoToBegin();!ItI_GSZM.IsAtEnd(); ++ItI_GSZM, ++ItM_GSZM )<br> {<br> if ( int(ItM_GSZM.Get()) ==A_label )<br> {<br> level=ItI_GSZM.Get();<br> volume_all[level]=volume_all[level]+1;<br> }<br> }<br><br> int mask_volume=0;<br> int level_mask_volume=0;<br><br> for(int k=minValue; k<(maxValue+1); k++){<br> if (mask_volume<volume_all[k])<br> {<br> mask_volume=volume_all[k];<br> level_mask_volume=k;<br> }<br> // std::cout << " levels : " << k <<" volume : "<<volume_all[k] << std::endl;<br> }<br><br> //std::cout << " volume mask : " << mask_volume <<" voxels"<< std::endl;<br> //std::cout << " level volume mask : " << level_mask_volume << std::endl;<br><br><br><br>//now we have an empty 2D-matrix of size (0,0). Refill it with one single command:<br> ZM.assign( num_of_col, vector<int>(num_of_row) ) ;<br><br>//std::cout << "MaskImage initialaise 1"<< std::endl;<br><br>int volume=0;<br>int Max_volume=0;<br>int Min_volume=mask_volume;<br>int N_object_connect=0;<br><br>for(int level_Grey=minValue; level_Grey<(maxValue+1); level_Grey++ )<br>{<br> //choose one levels of Grey (D levels)<br> thresholdFilterToBinary_GSZM->SetInput(image_GSZM);<br> thresholdFilterToBinary_GSZM->SetLowerThreshold(int(level_Grey));<br> thresholdFilterToBinary_GSZM->SetUpperThreshold(int(level_Grey));<br> thresholdFilterToBinary_GSZM->SetInsideValue(int(level_Grey));<br> thresholdFilterToBinary_GSZM->SetOutsideValue(0);<br> thresholdFilterToBinary_GSZM->Update();<br><br> //look at the connectivity between voxels each zone is save in a different label<br> connected_GSZM->SetInput(thresholdFilterToBinary_GSZM->GetOutput());<br> connected_GSZM->SetFullyConnected(1);<br> connected_GSZM->Update();<br><br> //std::cout << "Number of objects for D= "<<level_Grey <<": "<< connected_GSZM->GetObjectCount() << std::endl;<br> N_object_connect=0;<br> N_object_connect=connected_GSZM->GetObjectCount();<br><br> //measure the volume of each labels<br> labelGeometryImageFilter_GSZM->SetInput(connected_GSZM->GetOutput() );<br> labelGeometryImageFilter_GSZM->SetIntensityInput(image_GSZM);<br> labelGeometryImageFilter_GSZM->Update();<br><br> //store in the Thibault matrix<br> for(int N=1; N<int(N_object_connect+1);N++)<br> {<br> volume=labelGeometryImageFilter_GSZM->GetVolume(N);<br> ZM[level_Grey][volume]++;<br> // std::cout << "Grey value of: "<<level_Grey+1<< " Matrix filling volume: " << volume <<" voxels"<< std::endl;<br> if ( Max_volume<volume ) {Max_volume=volume;}<br> if ( Min_volume>volume ) {Min_volume=volume;}<br><br> }<br> //inputImage->DisconnectPipeline();<br>}<br><br>//std::cout << "GSZM create"<< std::endl;<br>std::cout << "max volume :"<< Max_volume<< std::endl;<br><br><br><br>//////////////////////////////////////////////////////////////////////////////////////////////////:<br><br>/////////////////////convert size zone matrix to an histogram//////<br><br>typedef itk::Statistics::SparseFrequencyContainer2 FrequencyContainerType_GSZM;<br>typedef FrequencyContainerType_GSZM::AbsoluteFrequencyType FrequencyType_GSZM;<br><br>typedef itk::Statistics::Histogram< float, FrequencyContainerType_GSZM> HistogramType_GSZM;<br><br> unsigned int numberOfComponents=2;<br> HistogramType_GSZM::Pointer histogram = HistogramType_GSZM::New();<br><br> HistogramType_GSZM::SizeType size(numberOfComponents);<br> size[0]=D;<br> size[1]=(Max_volume-Min_volume);<br><br> HistogramType_GSZM::MeasurementVectorType lowerBound(numberOfComponents);<br> lowerBound[0]=0;<br> lowerBound[1]=0;<br><br> HistogramType_GSZM::MeasurementVectorType upperBound(numberOfComponents);<br> upperBound[0]=(maxValue-minValue);<br> upperBound[1]=Max_volume-1;<br><br> histogram->SetMeasurementVectorSize(numberOfComponents);<br> histogram->Initialize(size, lowerBound, upperBound );<br> histogram->SetToZero();<br><br> // std::cout << "Histogramme created: middle "<< std::endl;<br><br> HistogramType_GSZM::IndexType index(numberOfComponents);<br> for(int level_Grey=minValue; level_Grey<(maxValue+1); level_Grey++ )<br> {<br> for(int N=0; N<(Max_volume);N++)<br> {<br> index[0]=(level_Grey-minValue); //the calulation add +1 to the two index<br> index[1]=N;<br> histogram->SetFrequencyOfIndex(index, static_cast<FrequencyType_GSZM>(ZM[level_Grey][N+1]));<br> /* if (N==(1) )<br> {<br> std::cout << "index: "<<index[0]<<" "<<index[1]<<" value: "<< Z[level_Grey][N]<< std::endl;<br> }*/<br> }<br> }<br> // std::cout << "Histogramme created: ok "<< std::endl;<br><br> //-------------------libére les vecteurs:<br><br> //ZM.clear();<br><br><br>//////////////////////////////////calcul of matrix size features////////////////<br>/////////////:based on the run lenght matrix////////////////:<br><br>typedef itk::Statistics::HistogramToRunLengthFeaturesFilter<HistogramType_GSZM> HistoSizeFilterType;<br><br>HistoSizeFilterType::Pointer hzf= HistoSizeFilterType::New();<br>hzf->SetInput(histogram);<br><br>//hzf->FastCalculationsOn();//just one direction<br><br>hzf->Update();<br><br>std::cout<<" GSZM: ok"<<std::endl;<br><br><br>//And then to acess the values<br></div><div><br> //GSZM:<br><br> fichier<<hzf->GetShortRunEmphasis()<<"\t";<br> fichier<<hzf->GetLongRunEmphasis()<<"\t";<br> fichier<<hzf->GetGreyLevelNonuniformity()<<"\t";<br> fichier<<hzf->GetRunLengthNonuniformity()<<"\t";<br> fichier<<hzf->GetLowGreyLevelRunEmphasis()<<"\t";<br> fichier<<hzf->GetHighGreyLevelRunEmphasis()<<"\t";<br> fichier<<hzf->GetShortRunLowGreyLevelEmphasis()<<"\t";<br> fichier<<hzf->GetShortRunHighGreyLevelEmphasis()<<"\t";<br> fichier<<hzf->GetLongRunLowGreyLevelEmphasis()<<"\t";<br> fichier<<hzf->GetLongRunHighGreyLevelEmphasis()<<"\t";<br><br></div><div>°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°</div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-04-21 19:31 GMT+02:00 Pati, Sarthak <span dir="ltr"><<a href="mailto:Sarthak.Pati@uphs.upenn.edu" target="_blank">Sarthak.Pati@uphs.upenn.edu</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div link="blue" vlink="purple" lang="EN-US">
<div class="m_7090690268976619214WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d">Hi Cyril,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d">Sure, that would be fine, as long as there is nothing missing.
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d">Thanks so much,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d">Sarthak
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Segoe UI",sans-serif;color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><b><span style="font-size:11.0pt;font-family:"Calibri",sans-serif">From:</span></b><span style="font-size:11.0pt;font-family:"Calibri",sans-serif"> Cyril Jaudet [mailto:<a href="mailto:drcjaudet@gmail.com" target="_blank">drcjaudet@gmail.com</a>]
<br>
<b>Sent:</b> 21/Apr/2017; Friday 13:31<br>
<b>To:</b> Pati, Sarthak <<a href="mailto:Sarthak.Pati@uphs.upenn.edu" target="_blank">Sarthak.Pati@uphs.upenn.edu</a>><br>
<b>Subject:</b> Re: Grey Level Size Zone Matrix<u></u><u></u></span></p><div><div class="h5">
<p class="MsoNormal"><u></u> <u></u></p>
<div>
<div>
<div>
<p class="MsoNormal">Hello Pati, <u></u><u></u></p>
</div>
<p class="MsoNormal">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.
<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">See you, <u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">Cyril Jaudet, PhD<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">Medical physicist<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">UZBrussel, Belgium<u></u><u></u></p>
</div>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
<div>
<p class="MsoNormal">2017-04-21 18:10 GMT+02:00 Pati, Sarthak <<a href="mailto:Sarthak.Pati@uphs.upenn.edu" target="_blank">Sarthak.Pati@uphs.upenn.edu</a>>:<u></u><u></u></p>
<blockquote style="border:none;border-left:solid #cccccc 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">Hi,</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif"> </span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">I saw your post on the ITK forums (<a href="http://public.kitware.com/pipermail/community/2015-January/008102.html" target="_blank">http://public.kitware.com/<wbr>pipermail/community/2015-<wbr>January/008102.html</a>)
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.</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif"> </span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">Thanks,</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">Sarthak Pati</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">Lead Software Developer</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">Center for Biomedical Image Computing & Analytics</span><u></u><u></u></p>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif">University of Pennsylvania</span><u></u><u></u></p>
<p class="MsoNormal"> <u></u><u></u></p>
</div>
</div>
</blockquote>
</div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
</div></div></div>
</div>
</blockquote></div><br></div>