[Insight-users] FFTRealToComplexConjugateImageFilter - bug +
fix + question
Kent Williams
kent at psychiatry.uiowa.edu
Wed Sep 20 10:18:09 EDT 2006
Jakub, thanks for your sleuthing. This -- or some version of this --
will go into ITK CVS.
I recently added a new function to the
itk::FFTComplexConjugateToRealFilter to get around the problem with
recovering the size of the real input from the size of the complex
conjugate:
void SetActualXDimensionIsOdd(bool isodd)
{
m_ActualXDimensionIsOdd = isodd;
}
The rationale is this: Previously, the actual X dimension of the real
image was attached as metadata to the output image of the Real to
Complex filter. The metadata would propogate down the pipeline from
input to output image, and the Complex to Real filter would use the
actual X dimension reported in metadata to size the outpt image. This
works fine as long as you're using a pipeline with a R->C ... C->R
structure, but one of the researchers here at Iowa was actually saving
the Complex image to disk, and image metadata doesn't get saved in the
image file.
This way, when the application reads in the complex image, the program
can tell the filter whether the X Dimension was odd, and the correctly
sized output image is used. I'd be the first to admit that this is a
little clumsy, but it's an inevitable property of FFT/IFFT. Actually
encountering this problem would be rare just because imaging
applications rarely use an odd number for the X size of images, but the
solution (more a workaround than a solution) is there now.
Jakub Bican wrote:
>
> 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
>
>
>------------------------------------------------------------------------
>
>/*=========================================================================
>
> 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
>
>
>------------------------------------------------------------------------
>
>/*=========================================================================
>
> 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
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>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