[Insight-developers] [Fwd: [Insight-users]
FFTRealToComplexConjugateImageFilter - bug + fix + question]
Kent Williams
kent at psychiatry.uiowa.edu
Wed Sep 20 10:03:17 EDT 2006
Jakub Bican posted this message about the itkFFT... filters on the
Insight Users list. I'm not an expert on ITK pipeline operation, but
what he says here seems plausible. Should I go ahead and check his
change into CVS ITK?
-------- Original Message --------
Subject: [Insight-users] FFTRealToComplexConjugateImageFilter - bug +
fix + question
Date: Sun, 17 Sep 2006 20:28:42 +0200
From: Jakub Bican <jakub.bican at matfyz.cz>
Reply-To: jakub.bican at matfyz.cz
To: Insight-users at itk.org
Hi,
after more than year, i am using FFT classes in ITK in my new project.
And i reached the same bug as one year ago:))
( see mail "FW: [Insight-users] FFTW library in ITK" from 21.5.2005 for
the bug description and testing code )
BUT TODAY - i found where is the problem... :))
The FFTRealToComplexConjugateImageFilter sets the LargestPossibleRegion
of the output, so that it reflets the halving effect of the FFT in the
first dimension (i.e. outputSize[0] = inputSize[0]/2+1). Then
ImageToImageFilter (or some other ancestor) propagates this size to the
RequestedRegion of the input image - and that is the problem, because
FFTWReal...Filter accesses buffer of the input image as if it is twice
as big.
So if the first upstream filter from the FFT filter is Update()d
manually, the BufferedRegion is equal to the LargestPossibleRegion and
everything works. But if the filter is updated due to upstream
propagation, the BufferedRegion is equal to RequestedRegion, that is
half of the LargestPossibleRegion. The FFT filter relies on that
BufferedRegion is equal to LargestPossibleRegion and reaches a serious
problem.
*THE SOLUTION*
add this to the FFTComplexConjugateFilter:
template < class TPixel , unsigned int Dimension >
void
FFTRealToComplexConjugateImageFilter < TPixel , Dimension >
::GenerateInputRequestedRegion()
{
Superclass::GenerateInputRequestedRegion();
ImageType::Pointer input = const_cast<ImageType *> (this->GetInput());
if ( !input )
{
return;
}
input->SetRequestedRegionToLargestPossibleRegion();
}
template < class TPixel , unsigned int Dimension >
void
FFTRealToComplexConjugateImageFilter < TPixel , Dimension >
::EnlargeOutputRequestedRegion(DataObject *output)
{
Superclass::EnlargeOutputRequestedRegion(output);
output->SetRequestedRegionToLargestPossibleRegion();
}
The second method is not absolutely neccessary in this case (so far i
know:).
The FFTWRealToComplexConjugateFilter should also check if the input data
size matches the assumed size.
*THE QUESTION*
Why this bug occurred only for some instances of the template? (Now, it
happens when the pixel type of the FFT filter is double.)
I am attaching the fixed files. I did not document it very much as the
whole class is not documented very much.
Regards,
Jakub
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkFFTRealToComplexConjugateImageFilter.h,v $
Language: C++
Date: $Date: 2004/11/01 16:29:56 $
Version: $Revision: 1.3 $
Copyright (c) 2002 Insight Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __itkFFTRealToComplexConjugateImageFilter_h
#define __itkFFTRealToComplexConjugateImageFilter_h
#include <itkImageToImageFilter.h>
#include <itkImage.h>
#include <complex>
namespace itk
{
/** /class FFTRealToComplexConjugateImageFilter
* /brief
*
* \ingroup
*/
template <class TPixel,unsigned int Dimension = 3>
class ITK_EXPORT FFTRealToComplexConjugateImageFilter :
public ImageToImageFilter< Image< TPixel , Dimension >,
Image< std::complex< TPixel > , Dimension > >
{
public:
/** Standard class typedefs.*/
typedef Image<TPixel,Dimension> TInputImageType;
typedef Image< std::complex< TPixel > , Dimension> TOutputImageType;
typedef FFTRealToComplexConjugateImageFilter Self;
typedef ImageToImageFilter< TInputImageType, TOutputImageType > Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> constPointer;
/** Run-time type information (and related methods). */
itkTypeMacro(FFTRealToComplexConjugateImageFilter, ImageToImageFilter);
/** Image type typedef support. */
typedef TInputImageType ImageType;
typedef typename ImageType::SizeType ImageSizeType;
virtual void GenerateOutputInformation(); // figure out allocation for output image
virtual void GenerateInputRequestedRegion();
virtual void EnlargeOutputRequestedRegion(DataObject *output);
virtual bool FullMatrix() = 0; // must be implemented in child
protected:
FFTRealToComplexConjugateImageFilter() {}
virtual ~FFTRealToComplexConjugateImageFilter(){}
private:
FFTRealToComplexConjugateImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFFTRealToComplexConjugateImageFilter.txx"
#endif
#endif
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkFFTRealToComplexConjugateImageFilter.txx,v $
Language: C++
Date: $Date: 2003/11/10 18:51:16 $
Version: $Revision: 1.2 $
Copyright (c) 2002 Insight Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __itkFFTRealToComplexConjugateImageFilter_txx
#define __itkFFTRealToComplexConjugateImageFilter_txx
#include "itkFFTRealToComplexConjugateImageFilter.h"
#include "itkMetaDataObject.h"
namespace itk
{
template < class TPixel , unsigned int Dimension >
void
FFTRealToComplexConjugateImageFilter < TPixel , Dimension >
::GenerateOutputInformation()
{
// call the superclass' implementation of this method
Superclass::GenerateOutputInformation();
//
// If this implementation returns a full result
// instead of a 'half-complex' matrix, then none of this
// is necessary
if(this->FullMatrix())
return;
// get pointers to the input and output
typename TInputImageType::ConstPointer inputPtr = this->GetInput();
typename TOutputImageType::Pointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
//
// This is all based on the same function in itk::ShrinkImageFilter
// ShrinkImageFilter also modifies the image spacing, but spacing
// has no meaning in the result of an FFT.
unsigned int i;
const typename TInputImageType::SizeType& inputSize
= inputPtr->GetLargestPossibleRegion().GetSize();
const typename TInputImageType::IndexType& inputStartIndex
= inputPtr->GetLargestPossibleRegion().GetIndex();
typename TOutputImageType::SizeType outputSize;
typename TOutputImageType::IndexType outputStartIndex;
//
// in 4.3.4 of the FFTW documentation, they indicate the size of
// of a real-to-complex FFT is N * N ... + (N /2+1)
// 1 2 d
// complex numbers.
// static_cast prob. not necessary but want to make sure integer
// division is used.
outputSize[0] = static_cast<unsigned int>(inputSize[0])/2 + 1;
outputStartIndex[0] = inputStartIndex[0];
for (i = 1; i < TOutputImageType::ImageDimension; i++)
{
outputSize[i] = inputSize[i];
outputStartIndex[i] = inputStartIndex[i];
}
//
// the halving of the input size hides the actual size of the input.
// to get the same size image out of the IFFT, need to send it as
// Metadata.
typedef typename TOutputImageType::SizeType::SizeValueType SizeScalarType;
itk::MetaDataDictionary &OutputDic=outputPtr->GetMetaDataDictionary();
itk::EncapsulateMetaData<SizeScalarType>(OutputDic,
std::string("FFT_Actual_RealImage_Size"),
inputSize[0]);
typename TOutputImageType::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize( outputSize );
outputLargestPossibleRegion.SetIndex( outputStartIndex );
outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
}
template < class TPixel , unsigned int Dimension >
void
FFTRealToComplexConjugateImageFilter < TPixel , Dimension >
::GenerateInputRequestedRegion()
{
// call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// get pointers to the inputs
ImageType::Pointer input = const_cast<ImageType *> (this->GetInput());
if ( !input )
{
return;
}
input->SetRequestedRegionToLargestPossibleRegion();
}
template < class TPixel , unsigned int Dimension >
void
FFTRealToComplexConjugateImageFilter < TPixel , Dimension >
::EnlargeOutputRequestedRegion(DataObject *output)
{
Superclass::EnlargeOutputRequestedRegion(output);
output->SetRequestedRegionToLargestPossibleRegion();
}
}
#endif
-------------- next part --------------
_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-developers
mailing list