[Insight-users] EqalizationFilter: Mini-Pipeline Problem
Karthik Krishnan
Karthik.Krishnan at kitware.com
Fri Aug 5 13:20:44 EDT 2005
Ok.. Looking at the code, I can see that you intend reusing the output
of the rescale filter. To correct myself, I just tried your code and it
works fine. The scalarToArrayCast filter does work on the rescaled
output. Are you having a problem using it ?
thanks
karthik
On Fri, 2005-08-05 at 12:58 -0400, Karthik Krishnan wrote:
> You should not have
>
> rescalePtr->GraftOutput( this->GetOutput() );
>
> If you do, the itInput and itOutput are going to point to the same image
> container data.
>
>
> hth
> karthik
>
> Torsten Schröer wrote:
>
> >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.
> >
> >
> >
> >
> > ------------------------------------------------------------------------
> >
> >------------------------------------------------------------------------
> >
> >// 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
> >
> >
> >------------------------------------------------------------------------
> >
> >// 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
> >
> >
> >------------------------------------------------------------------------
> >
> >_______________________________________________
> >Insight-users mailing list
> >Insight-users at itk.org
> >http://www.itk.org/mailman/listinfo/insight-users
> >
> >
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list