[Insight-users] Filter running out of memory

[Bradley Lowekamp]

Hello,

I am not entirely sure what your filter is trying to do, nor what's the problem you are encountering.

I have an itk module which has a filter to compute the GLCM for a neighborhood around each pixel. This may be what you are looking for, or at least give you some ideas on how to use the components.

As it has been a while since I last used this filter, and there is a lack of rigorous testing, I'd recommend writing a test to verify the results for your usage before relaying on the output.

Good luck,
Brad

[1] https://github.com/blowekamp/itkTextureAnalysis/blob/master/include/itkTextureFeatureImageFilter.h

On Aug 19, 2013, at 8:46 PM, Fabian Torres <dae.wong at gmail.com> wrote:

> Hi all. I am trying to implement a filter that gives local texture analysis. For this I implement a itkImagetoImageFilter that calculates coocurrence texture descriptor for each pixel using itkScalarImageToCoocurrenceMatrizFilter.
> 
> I know this is a slow implementation, but the real problem I have is that when I try to use this with a 3D image, even a small one of 43x35x61, the computer runs out of memory. I do not know if I my implementation is wrong. But it seems that for each pixel I´m using a lot of memory.
> 
> I hope someone can help me with this. Here I leave the code for the GenerateData() function.
> 
> Thanks.
> 
> this->AllocateOutputs();
> 
>  ProgressReporter progress( this, 0, 
>        this->GetInput()->GetRequestedRegion().GetNumberOfPixels(), 100 );
> 
>  typedef MirrorPadImageFilter<InputImageType,InputImageType> PadFilterType;
>  typename PadFilterType::Pointer padFilter = PadFilterType::New();
> 
>  int boundSize = floor(this->m_RegionSize[1]/2);
> 
>  RegionSizeType bound;
>  bound.Fill(boundSize);
> 
>  padFilter->SetPadBound(bound);
>  padFilter->SetInput(this->GetInput());
>  padFilter->Update();
> 
>  typename InputImageType::Pointer padImage = padFilter->GetOutput();
> 
>  typedef Statistics::ScalarImageToCooccurrenceMatrixFilter<InputImageType>
>          GLCMGeneratorType;
>  typename GLCMGeneratorType::Pointer glcmGenerator = GLCMGeneratorType::New();
> 
>  glcmGenerator->SetNumberOfBinsPerAxis(this->m_NumberofBins);
>  glcmGenerator->SetPixelValueMinMax(this->m_PixelMinValue,m_PixelMaxValue);
> 
>  typedef typename InputImageType::OffsetType OffsetType;
> 
>  OffsetType offset1;
>  offset1[0] = 1;
>  offset1[1] = 0;
> 
>  OffsetType offset2;
>  offset2[0] = 1;
>  offset2[1] = -1;
> 
>  OffsetType offset3;
>  offset3[0] = 0;
>  offset3[1] = -1;
> 
>  OffsetType offset4;
>  offset4[0] = -1;
>  offset4[1] = -1;
> 
>  typedef typename GLCMGeneratorType::OffsetVector OffsetVectorType;
>  typename OffsetVectorType::Pointer offsets = OffsetVectorType::New();  
> 
>  offsets->reserve(4);
>  offsets->InsertElement(0,offset1);
>  offsets->InsertElement(1,offset2);
>  offsets->InsertElement(2,offset3);
>  offsets->InsertElement(3,offset4);
> 
>  glcmGenerator->SetOffsets(offsets);
> 
>  typedef typename GLCMGeneratorType::HistogramType HistogramType;
>  typedef Statistics::HistogramToTextureFeaturesFilter<HistogramType> 
>          TextureFeaturesType;
>  typename TextureFeaturesType::Pointer textureFeatures;
> 
>  if(this->m_TextureFeature<=8 && this->m_TextureFeature>=0){  
>      textureFeatures = TextureFeaturesType::New();
>  }else{
>      textureFeatures = NULL;
>  }
> 
>  typedef ExtractImageFilter<InputImageType,InputImageType> ExtractFilterType;
>  typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
> 
>  typename InputImageType::RegionType region;
>  typename InputImageType::IndexType regionIndex;
> 
>  region.SetSize(this->m_RegionSize);
> 
>  extractFilter->SetInput(padImage);
> 
>  ImageRegionIteratorType itPad(padImage, padImage->GetLargestPossibleRegion());
>  itPad.GoToBegin();
> 
>  ImageRegionIteratorType itOut(this->GetOutput(), 
>          this->GetOutput()->GetLargestPossibleRegion());
>  itOut.GoToBegin();
> 
>  while(!itPad.IsAtEnd()){
> 
>        regionIndex = itPad.GetIndex();
>        region.SetIndex(regionIndex);
> 
>        if(padImage->GetLargestPossibleRegion().IsInside(region)){ 
> 
>            extractFilter->SetExtractionRegion(region);
>            extractFilter->UpdateLargestPossibleRegion();
> 
>            glcmGenerator->SetInput(extractFilter->GetOutput());
>            glcmGenerator->UpdateLargestPossibleRegion();
> 
>            if(this->m_TextureFeature == 0){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::Energy));
>            }else if(this->m_TextureFeature == 1){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::Entropy));
>            }else if(this->m_TextureFeature == 2){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::Correlation));
>            }else if(this->m_TextureFeature == 3){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::InverseDifferenceMoment));
>            }else if(this->m_TextureFeature == 4){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::Inertia)); 
>            }else if(this->m_TextureFeature == 5){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::ClusterShade));
>            }else if(this->m_TextureFeature == 6){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::ClusterProminence));
>            }else if(this->m_TextureFeature == 7){
>                textureFeatures->SetInput(glcmGenerator->GetOutput());
>                textureFeatures->UpdateLargestPossibleRegion();
>                this->GetOutput()->SetPixel(itOut.GetIndex(),
>                        textureFeatures->GetFeature(TextureFeaturesType::HaralickCorrelation));
>            }else if(this->m_TextureFeature == 8){
> 
>                typename HistogramType::ConstPointer hist = glcmGenerator->GetOutput();
> 
>                typename HistogramType::ConstIterator it = hist->Begin();
> 
>                float mean = 0;
>                while(it != hist->End()){
>                    mean += it.GetFrequency();
>                    ++it;
>                }
> 
>                float nBins = hist->GetSize(0);
>                nBins *= nBins;
> 
>                mean /= nBins;
> 
>                float variance = 0;
>                it = hist->Begin();
>                while( it != hist->End()){
>                    variance += pow(it.GetFrequency()-mean,2);
>                    ++it;
>                }
> 
>                variance /= (nBins - 1);
> 
>                this->GetOutput()->SetPixel(itOut.GetIndex(),variance);
> 
>                hist = NULL;
> 
>            }
> 
>            ++itOut;
> 
>            progress.CompletedPixel();
> 
>        }
> 
>        ++itPad; 
> _____________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
> 
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
> 
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users