[Insight-users] FFTRealToComplexConjugateImageFilter - bug + fix
+ question
Jakub Bican
jakub.bican at matfyz.cz
Thu Sep 21 15:04:58 EDT 2006
Hi Kent,
i would vote for preserving both ways - MetaData + SetActual...() method. The
second one is good in situation you described, while the first one is suitable
for example in case, when some filter in the pipeline, that is between forward
and inverse FFT filters changes the size of data - for example crops some
higher frequencies and changes spacing to subsample the image. If there is
information in metadata, it may be changed by the filter as well. If there is
only the new method, noone would know if the cropped image has to have odd or
even x size.
One more way would be to enable a switch, that would force FFT filters to
FullMatrix mode. Then nothing is lost and nothing extra has to be passed (but
more memory and time is consumed to store and process the redundant data).
I feel a lot of controversy in all of these approaches. The problem is that
the image is transformed from one space (euler basis) to other space (fourier
basis), while the image properties like size, origin and spacing still relate
to the original space (they "speak in terms of original space"). This is the
difference between FFT filters and other filters.
{{forgive me the terms please - i am not sure with english terms and i do not
have the dictionary here now}}
Some general approach would be suggested, to overcome the similar (or even
worse) problems in case of transformations to other spaces and bases like
wavelets, moments, etc.
SO, lets say, what is the attitude of the community to this problem. Am i
wrong? Is it unimportant or too difficult to implement in current itk so that
it would be always work-arounded(huh:(), like the FFT is now?
uff, i've got enough today. regards,
Jakub
PS: Some by-the-way question:
What is the state of automatic passing of metadata through the pipeline now?
---------- Original Message -----------
From: Kent Williams <kent at psychiatry.uiowa.edu>
To: jakub.bican at matfyz.cz
Cc: Insight-users at itk.org
Sent: Wed, 20 Sep 2006 09:18:09 -0500
Subject: Re: [Insight-users] FFTRealToComplexConjugateImageFilter - bug +
fix + question
> 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
> >
> >
------- End of Original Message -------
More information about the Insight-users
mailing list