[Insight-users] FFTRealToComplexConjugateImageFilter - bug +
fix + question
Kent Williams
kent at psychiatry.uiowa.edu
Thu Sep 28 19:22:00 EDT 2006
The FullMatrix method is informational. It's meant to be redefined in
subclasses -- i.e. for VNL, FFTW etc -- so that you can determine
whether the implementation returns a full matrix, or a complex
conjugate, which is a half matrix. Some FFT implementations might allow
both full-matrix or complex conjugate operations, but I don't know about
them. At this point, ITK CVS does still set the MetaDataDictionary
variable with the actual first dimension, but I'd keep track of it at
the application level and call the ActualXDimensionIsOdd method to clue
in the ComplexConjugateToReal filter.
I actually don't know anything more about using FFT/IFFT than I needed
to to write the filters. If you think resampling the complex result of
FFT, then using the IFFT filter generates meaningful results, you
already know more about it than I do.
I think if you want to understand the filters, you should read the
source -- what I did was pretty straightforward, and it's a good thing
if you're using ITK to understand how ITK filters work. If you have
ideas about how to improve the filters, you could try implementing them
and submit them to the Insight Journal for inclusion in ITK.
http://www.insight-journal.org/
Jakub Bican wrote:
> 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 -------
>
> _______________________________________________
> 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