[Insight-users] EqalizationFilter: Mini-Pipeline Problem

[Torsten Schröer]

Hi everybody!

Please take a look at the attached source of EqualizationFilter. I try
to set up a mini Pipeline inside this Filter connecting
IntensityRescaleFilter's output to ScalarToArrayCastFilter's input.
Something must go wrong there:

    // Rescaling
    RescaleIntensityFilterType::Pointer rescalePtr =
      RescaleIntensityFilterType::New();
    rescalePtr->GraftOutput( this->GetOutput() );
    rescalePtr->SetInput( inputPtr );
    rescalePtr->SetOutputMinimum(InputPixelTraits::Zero);
    rescalePtr->SetOutputMaximum(InputPixelTraits::max());
    rescalePtr->Update();

    // Convert Image
    ScalarToArrayCastFilterType::Pointer convertPtr =
      ScalarToArrayCastFilterType::New();
    convertPtr->SetInput( rescalePtr->GetOutput() );
    convertPtr->Update();

It seems as ScalarToArrayCastFilter's Input refers to the global Input
instead of the RescaleFilter's output. I tried to find more information
on excecuting mini-pipelines and looked at OtsuMultipleThresholds and
similar filters that use an mini pipeline whitout success.

This filter should work on unsigned data images but not jet on signed
data (work in progress).

When I put the Rescale Filter outside the mini-pipeline, just before the
EqualizationFilter, everything works properly (See attached histo).

If You think the filter might be useful, feel free to use it, but do not
forget about rescaling before executing the filter ;)

Any ideas appreciated,
Thanks,

Torsten.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: histo.jpg
Type: image/jpeg
Size: 9700 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20050804/8d2581df/histo.jpg
-------------- next part --------------
// itkEqualizationImageFilter.h
// Torsten Schröer
// Created: 04.08.2005
//
//////////////////////////////////////////////////////////////////////

#ifndef __ITK_EQUALIZATION_IMAGE_FILTER_H
#define __ITK_EQUALIZATION_IMAGE_FILTER_H

#include "itkImage.h"
#include "itkFixedArray.h"
#include "itkHistogram.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "itkImageToHistogramGenerator.h"
#include "itkImageToImageFilter.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "itkProgressAccumulator.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkScalarToArrayCastImageFilter.h"

namespace itk
{

/** 
 * \class EqualizationImageFilter performs a (very) simple histogram equalization.
 */
template <class TInputImage>
class ITK_EXPORT EqualizationImageFilter : 
  public ImageToImageFilter<TInputImage, TInputImage>
{
public:
  /** Standard class typedefs. */
  typedef TInputImage InputImageType;
  typedef TInputImage OutputImageType;

  typedef EqualizationImageFilter                       Self;
  typedef ImageToImageFilter<TInputImage, TInputImage>  Superclass;
  typedef SmartPointer<Self>                            Pointer;
  typedef SmartPointer<const Self>                      ConstPointer;
  
  /* Method for creation through the object factory. */
  itkNewMacro(Self);  

  const char* GetClassName () { return "EqualizationImageFilter"; } 

  /** Typedefs regarding input image */
  typedef typename InputImageType::ConstPointer         InputImagePointer;
  typedef typename InputImageType::PixelType            InputPixelType;
  typedef itk::NumericTraits<InputPixelType>            InputPixelTraits;

  /** Typedefs regarding output image */
  typedef typename OutputImageType::Pointer             OutputImagePointer;
  typedef typename OutputImageType::PixelType           OutputPixelType;
  typedef itk::NumericTraits<OutputPixelType>           OutputPixelTraits;

  /** Typedefs regarding image iterators */
  typedef itk::ImageRegionConstIterator<InputImageType> InputIteratorType;
  typedef itk::ImageRegionIterator<OutputImageType>     OutputIteratorType;

  /** Typedefs for HistogramGenerator */
  typedef itk::FixedArray< InputPixelType, 1 >                              MeasurementVectorType;
  typedef itk::Image< MeasurementVectorType, 3 >                            ArrayImageType;
  typedef itk::ScalarToArrayCastImageFilter<InputImageType, ArrayImageType> ScalarToArrayCastFilterType;
  typedef itk::Statistics::ImageToHistogramGenerator<ArrayImageType>        HistogramGeneratorType;
  typedef HistogramGeneratorType::HistogramType                             HistogramType;
  typedef typename HistogramType::SizeType                                  HistogramSizeType;
  
  /** Typedefs regarding rescaling */
  typedef itk::RescaleIntensityImageFilter<TInputImage, TInputImage> RescaleIntensityFilterType;

protected:
  EqualizationImageFilter();
  ~EqualizationImageFilter() {};

  // Print Self
  void PrintSelf(std::ostream& os, Indent indent) const;

	// Override since the filter needs all the data for the algorithm
	void GenerateInputRequestedRegion();

	// Override since the filter produces the entire dataset
	void EnlargeOutputRequestedRegion(DataObject *output);

  // The algorithm itself.
	void GenerateData();

private:
  EqualizationImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&);     //purposely not implemented

  HistogramGeneratorType::Pointer m_Generator;
  const HistogramType*            m_Histogram;
  HistogramType::Pointer          m_OutputHistogram;
};

  
} // end namespace itk
  
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkEqualizationImageFilter.txx"
#endif
  
#endif //__ITK_EQUALIZATION_IMAGE_FILTER_H
-------------- next part --------------
// itkEqualizationImageFilter.txx
// Torsten Schröer
// Created: 04.08.2005
//
//////////////////////////////////////////////////////////////////////


#ifndef _ITK_EQUALIZATION_IMAGE_FILTER_TXX
#define _ITK_EQUALIZATION_IMAGE_FILTER_TXX

#include "itkEqualizationImageFilter.h"

namespace itk
{

  template <class TInputImage>
  EqualizationImageFilter<TInputImage>
  ::EqualizationImageFilter()
  {
	  m_Generator       = HistogramGeneratorType::New();
    
    m_OutputHistogram = HistogramType::New();
  }

  template <class TInputImage>
  void 
  EqualizationImageFilter<TInputImage>
  ::PrintSelf( std::ostream& os, Indent indent ) const
  {
    Superclass::PrintSelf(os,indent);
    os << indent << "Equalization Image Filter."
	   << std::endl;
  }

  template <class TInputImage>
  void
  EqualizationImageFilter<TInputImage>
  ::GenerateInputRequestedRegion()
  {
	  Superclass::GenerateInputRequestedRegion();
	  if ( this->GetInput() ) {
		  InputImageType* image =
			  const_cast< InputImageType* >( this->GetInput() );
		  image->SetRequestedRegionToLargestPossibleRegion();
	  }
  }

  template <class TInputImage>
  void
  EqualizationImageFilter<TInputImage>
  ::EnlargeOutputRequestedRegion(DataObject *output)
  {
	  Superclass::EnlargeOutputRequestedRegion(output);
	  output->SetRequestedRegionToLargestPossibleRegion();
  }

  template <class TInputImage>
  void 
  EqualizationImageFilter<TInputImage>
  ::GenerateData()
  {
    unsigned int i = 0;

    // Get the image input and output pointers
	  InputImagePointer  inputPtr  = this->GetInput();
	  OutputImagePointer outputPtr = this->GetOutput();

    // Zero the output
	  outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
	  outputPtr->Allocate();
	  outputPtr->FillBuffer( OutputPixelTraits::Zero );
    
    // Compute images Size for scaling factor
    InputImageType::SizeType sz = inputPtr->GetRequestedRegion().GetSize();
    unsigned int imageSize = 1;
    for (i=0; i<inputPtr->GetImageDimension(); i++) {
      imageSize *= static_cast<unsigned int>(sz[i]);
    }

    // Prepare Progress Reporter
    ProgressReporter filterProgress(this, 0, imageSize);

    // Rescaling to whole range
    RescaleIntensityFilterType::Pointer rescalePtr = RescaleIntensityFilterType::New();
    rescalePtr->GraftOutput( this->GetOutput() );
    rescalePtr->SetInput( inputPtr );
    rescalePtr->SetOutputMinimum(InputPixelTraits::Zero);
    rescalePtr->SetOutputMaximum(InputPixelTraits::max());
    rescalePtr->Update();

    /**
    itk::ImageFileWriter<InputImageType>::Pointer writer = itk::ImageFileWriter<InputImageType>::New();
    writer->SetFileName("test.mhd");
    writer->SetInput(rescalePtr->GetOutput());
    writer->Update();
    */

    // Convert Image
    ScalarToArrayCastFilterType::Pointer convertPtr = ScalarToArrayCastFilterType::New();
    convertPtr->SetInput( rescalePtr->GetOutput() );
    convertPtr->Update();

    // Compute Histogram
    m_Generator->SetInput( convertPtr->GetOutput() );
    HistogramType::SizeType histoSize = {InputPixelTraits::max()};
    m_Generator->SetNumberOfBins(histoSize);
    m_Generator->SetMarginalScale(10.0);
	  m_Generator->Compute();

    m_Histogram = m_Generator->GetOutput();

    m_OutputHistogram->Initialize(histoSize);
    m_OutputHistogram->SetFrequency( 0, 0 /*m_Histogram->GetFrequency(0)*/ );
    // ofstream output("Histogram.dat");
    for (i=1; i<m_Histogram->Size(); i++) {
      HistogramType::FrequencyType fr0 = m_OutputHistogram->GetFrequency(i-1);
      HistogramType::FrequencyType fr1 = m_Histogram->GetFrequency(i);
      m_OutputHistogram->SetFrequency( i, (fr0 + fr1) );
      // output << i << "\t" << (fr0+fr1) << std::endl;
    }
    // output.close();

    InputIteratorType  itInput(inputPtr, inputPtr->GetRequestedRegion());
    OutputIteratorType itOutput(outputPtr, outputPtr->GetRequestedRegion());

    itInput.GoToBegin();
    itOutput.GoToBegin();

    while (!itInput.IsAtEnd() && !itOutput.IsAtEnd()) {
      OutputPixelType op = static_cast<OutputPixelType>(
        (m_Histogram->Size() - 1) * m_OutputHistogram->GetFrequency(itInput.Get()) / imageSize );
      itOutput.Set(op);
      ++itInput;
      ++itOutput;
      filterProgress.CompletedPixel();
    }

    this->GraftOutput( this->GetOutput() );
  }

} // end namespace itk

#endif