[Insight-users] FFTRealToComplexConjugateImageFilter - bug + fix + question

[Jakub Bican]

Skipped content of type multipart/alternative-------------- 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