[Insight-users] EqalizationFilter: Mini-Pipeline Problem
Torsten Schröer
torsten.schroeer at igd.fraunhofer.de
Sat Aug 6 06:07:00 EDT 2005
Hello Karthik,
thanks for testing and help! I just realized, that Iterator "itInput"
pointed to the global input instead of rescaling filters output. Argh ...
Now I fixed that (see attached code) and everything works as expected.
The problem I had occured, because I loaded an image as unsigned short,
which used only grayvalue range 0 to 450. Then rescaling, computation of
equalization and applying the result to global input leads to uncorrect
output data.
This filter fits my needs now, but still supports only <unsigned char>
and <unsigned short> Images. with other types m_Generator->Compute() fails.
Thanks again for your help,
Torsten.
Karthik Krishnan schrieb:
> 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
>
>
-------------- 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()
{
int i = 0;
InputPixelType minimum = InputPixelTraits::min();
InputPixelType maximum = InputPixelTraits::max();
InputPixelType range = maximum - minimum;
// 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->SetInput( inputPtr );
rescalePtr->SetOutputMinimum( minimum );
rescalePtr->SetOutputMaximum( maximum );
rescalePtr->Update();
// Convert Image
ScalarToArrayCastFilterType::Pointer convertPtr = ScalarToArrayCastFilterType::New();
convertPtr->SetInput( rescalePtr->GetOutput() );
convertPtr->Update();
// Compute Histogram
m_Generator->SetInput( convertPtr->GetOutput() );
HistogramType::SizeType histoSize = { range };
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, minimum );
for (i=1; i<static_cast<int>(range); i++) {
HistogramType::FrequencyType fr0 = m_OutputHistogram->GetFrequency( i-1 );
HistogramType::FrequencyType fr1 = m_Histogram->GetFrequency( i );
m_OutputHistogram->SetFrequency( i, (fr0 + fr1) );
}
InputIteratorType itInput( rescalePtr->GetOutput(), rescalePtr->GetOutput()->GetRequestedRegion() );
OutputIteratorType itOutput( outputPtr, outputPtr->GetRequestedRegion() );
itInput.GoToBegin();
itOutput.GoToBegin();
while (!itInput.IsAtEnd() && !itOutput.IsAtEnd()) {
InputPixelType inp = itInput.Get();
OutputPixelType op = static_cast<OutputPixelType>(
(m_Histogram->Size() - 1) * m_OutputHistogram->GetFrequency( minimum + inp ) / imageSize );
itOutput.Set( op );
++itInput;
++itOutput;
filterProgress.CompletedPixel();
}
this->GraftOutput( this->GetOutput() );
}
} // end namespace itk
#endif
More information about the Insight-users
mailing list