[Insight-users] EqalizationFilter: Mini-Pipeline Problem
Karthik Krishnan
Karthik.Krishnan at kitware.com
Fri Aug 5 12:58:17 EDT 2005
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
>
>
More information about the Insight-users
mailing list