[Insight-users] Filter running out of memory

[Fabian Torres]

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;