Commit f80b3799 authored by Julien Michel's avatar Julien Michel

COMP: Removing unused itk fork filter in PendingPatches, and test using it

parent 356050ce
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkFFTComplexToComplexImageFilter_h
#define __itkFFTComplexToComplexImageFilter_h
#include "itkImageToImageFilter.h"
#include <complex>
namespace itk
{
/** \class FFTComplexToComplexImageFilter
*
* \brief Implements an API to enable the Fourier transform or the inverse
* Fourier transform of images with complex valued voxels to be computed.
*
* \ingroup FourierTransform
*
* \author Simon K. Warfield simon.warfield@childrens.harvard.edu
*
* \note Attribution Notice. This research work was made possible by
* Grant Number R01 RR021885 (PI Simon K. Warfield, Ph.D.) from
* the National Center for Research Resources (NCRR), a component of the
* National Institutes of Health (NIH). Its contents are solely the
* responsibility of the authors and do not necessarily represent the
* official view of NCRR or NIH.
*
* This class was taken from the Insight Journal paper:
* http://hdl.handle.net/1926/326
*
* \ingroup FourierTransform
*
* \sa ForwardFFTImageFilter
* \ingroup ITKReview
*
* \ingroup OTBITKPendingPatches
*/
template< class TImage >
class FFTComplexToComplexImageFilter:
public ImageToImageFilter< TImage, TImage >
{
public:
/** Input and output image types. */
typedef TImage ImageType;
typedef TImage InputImageType;
typedef TImage OutputImageType;
/** Standard class typedefs. */
typedef FFTComplexToComplexImageFilter Self;
typedef ImageToImageFilter< InputImageType, OutputImageType > Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
itkStaticConstMacro(ImageDimension, unsigned int,
InputImageType::ImageDimension);
/** Run-time type information (and related methods). */
itkTypeMacro(FFTComplexToComplexImageFilter, ImageToImageFilter);
/** Customized object creation methods that support configuration-based
* selection of FFT implementation.
*
* Default implementation is FFTW.
*/
static Pointer New(void);
/** Transform Direction */
typedef enum {
DIRECT = 1,
INVERSE
} TransformDirectionType;
/** Image type typedef support. */
typedef typename ImageType::SizeType ImageSizeType;
/** Set/Get the direction in which the transform will be applied.
* By selecting DIRECT, this filter will perform a direct Fourier Transform,
* By selecting INVERSE, this filter will perform an inverse Fourier Transform,
*/
itkSetMacro(TransformDirection, TransformDirectionType);
itkGetConstMacro(TransformDirection, TransformDirectionType);
protected:
FFTComplexToComplexImageFilter() {}
virtual ~FFTComplexToComplexImageFilter(){}
/** methods needed for the image filter pipeline */
virtual void GenerateOutputInformation(); // figure out allocation for output
// image
virtual void GenerateInputRequestedRegion();
virtual bool FullMatrix() = 0; // must be implemented in child
private:
FFTComplexToComplexImageFilter(const Self &); //purposely not implemented
void operator=(const Self &); //purposely not implemented
TransformDirectionType m_TransformDirection;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFFTComplexToComplexImageFilter.hxx"
#endif
#endif
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
/**
*
* Attribution Notice. This research work was made possible by
* Grant Number R01 RR021885 (PI Simon K. Warfield, Ph.D.) from
* the National Center for Research Resources (NCRR), a component of the
* National Institutes of Health (NIH). Its contents are solely the
* responsibility of the authors and do not necessarily represent the
* official view of NCRR or NIH.
*
* This class was taken from the Insight Journal paper:
* http://hdl.handle.net/1926/326
*
*/
#ifndef __itkFFTComplexToComplexImageFilter_hxx
#define __itkFFTComplexToComplexImageFilter_hxx
#include "itkFFTComplexToComplexImageFilter.h"
#if defined( ITK_USE_FFTWD ) || defined( ITK_USE_FFTWF )
#include "itkFFTWComplexToComplexImageFilter.h"
#endif
namespace itk
{
#if defined( ITK_USE_FFTWD ) || defined( ITK_USE_FFTWF )
template< class TImage >
class FFTWComplexToComplexImageFilter;
#endif
template< class TImage >
typename FFTComplexToComplexImageFilter< TImage >::Pointer
FFTComplexToComplexImageFilter< TImage >
::New(void)
{
Pointer smartPtr = ::itk::ObjectFactory< Self >::Create();
#ifdef ITK_USE_FFTWD
if ( smartPtr.IsNull() )
{
if ( typeid( typename ImageType::PixelType::value_type ) == typeid( double ) )
{
smartPtr = dynamic_cast< Self * >(
FFTWComplexToComplexImageFilter< TImage >::New().GetPointer() );
}
}
#endif
#ifdef ITK_USE_FFTWF
if ( smartPtr.IsNull() )
{
if ( typeid( typename ImageType::PixelType::value_type ) == typeid( float ) )
{
smartPtr = dynamic_cast< Self * >(
FFTWComplexToComplexImageFilter< TImage >::New().GetPointer() );
}
}
#endif
return smartPtr;
}
template< class TImage >
void
FFTComplexToComplexImageFilter< TImage >::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 InputImageType::ConstPointer inputPtr = this->GetInput();
typename OutputImageType::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. For an IFFT, since the
// spacing is propagated to the complex result, we can use the spacing
// from the input to propagate back to the output.
unsigned int i;
const typename InputImageType::SizeType & inputSize =
inputPtr->GetLargestPossibleRegion().GetSize();
const typename InputImageType::IndexType & inputStartIndex =
inputPtr->GetLargestPossibleRegion().GetIndex();
typename OutputImageType::SizeType outputSize;
typename OutputImageType::IndexType outputStartIndex;
//
// Size of output FFT:C2C is the same as input
//
outputSize[0] = inputSize[0];
outputStartIndex[0] = inputStartIndex[0];
for ( i = 1; i < OutputImageType::ImageDimension; i++ )
{
outputSize[i] = inputSize[i];
outputStartIndex[i] = inputStartIndex[i];
}
typename OutputImageType::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize(outputSize);
outputLargestPossibleRegion.SetIndex(outputStartIndex);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
}
template< class TImage >
void
FFTComplexToComplexImageFilter< TImage >::GenerateInputRequestedRegion()
{
Superclass::GenerateInputRequestedRegion();
// get pointers to the input and output
typename InputImageType::Pointer inputPtr =
const_cast< InputImageType * >( this->GetInput() );
inputPtr->SetRequestedRegionToLargestPossibleRegion();
}
}
#endif
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkFFTWComplexToComplexImageFilter_h
#define __itkFFTWComplexToComplexImageFilter_h
#include "itkFFTComplexToComplexImageFilter.h"
#include "itkFFTWCommon.h"
namespace itk
{
/** \class FFTWComplexToComplexImageFilter
* \brief Implements an API to enable the Fourier transform or the inverse
* Fourier transform of images with complex valued voxels to be computed using
* either FFTW from MIT or the FFTW interface in Intel MKL.
* This filter is multithreaded and supports input images with sizes which are not
* a power of two.
*
* This code was contributed in the Insight Journal paper:
* "FFT Complex to Complex filters and helper classes"
* by Warfield S.
* http://hdl.handle.net/1926/326
* http://www.insight-journal.org/browse/publication/128
*
* \author Simon K. Warfield simon.warfield@childrens.harvard.edu
*
* \note Attribution Notice. This research work was made possible by
* Grant Number R01 RR021885 (PI Simon K. Warfield, Ph.D.) from
* the National Center for Research Resources (NCRR), a component of the
* National Institutes of Health (NIH). Its contents are solely the
* responsibility of the authors and do not necessarily represent the
* official view of NCRR or NIH.
*
* \ingroup FourierTransform
* \ingroup MultiThreaded
* \ingroup ITKReview
*
* \sa FFTWGlobalConfiguration
*
* \ingroup OTBITKPendingPatches
*/
template< class TImage >
class ITK_EXPORT FFTWComplexToComplexImageFilter:
public FFTComplexToComplexImageFilter< TImage >
{
public:
typedef FFTWComplexToComplexImageFilter Self;
typedef FFTComplexToComplexImageFilter< TImage > Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
/** Standard class typedefs. */
typedef TImage ImageType;
typedef typename ImageType::PixelType PixelType;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
//
// the proxy type is a wrapper for the fftw API
// since the proxy is only defined over double and float,
// trying to use any other pixel type is inoperative, as
// is trying to use double if only the float FFTW version is
// configured in, or float if only double is configured.
//
typedef typename fftw::Proxy< typename PixelType::value_type > FFTWProxyType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(FFTWComplexToComplexImageFilter,
FFTComplexToComplexImageFilter);
itkStaticConstMacro(ImageDimension, unsigned int,
ImageType::ImageDimension);
/** Image type typedef support. */
typedef typename ImageType::SizeType ImageSizeType;
//
// these should be defined in every FFT filter class
virtual bool FullMatrix();
/**
* Set/Get the behavior of wisdom plan creation. The default is
* provided by FFTWGlobalConfiguration::GetPlanRigor().
*
* The parameter is one of the FFTW planner rigor flags FFTW_ESTIMATE, FFTW_MEASURE,
* FFTW_PATIENT, FFTW_EXHAUSTIVE provided by FFTWGlobalConfiguration.
* /sa FFTWGlobalConfiguration
*/
virtual void SetPlanRigor( const int & value )
{
// use that method to check the value
FFTWGlobalConfiguration::GetPlanRigorName( value );
if( m_PlanRigor != value )
{
m_PlanRigor = value;
this->Modified();
}
}
itkGetConstReferenceMacro( PlanRigor, int );
void SetPlanRigor( const std::string & name )
{
this->SetPlanRigor( FFTWGlobalConfiguration::GetPlanRigorValue( name ) );
}
protected:
FFTWComplexToComplexImageFilter()
{
m_PlanRigor = FFTWGlobalConfiguration::GetPlanRigor();
}
virtual ~FFTWComplexToComplexImageFilter()
{
}
virtual void UpdateOutputData(DataObject *output);
virtual void BeforeThreadedGenerateData();
void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId );
void PrintSelf(std::ostream & os, Indent indent) const;
private:
FFTWComplexToComplexImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
bool m_CanUseDestructiveAlgorithm;
int m_PlanRigor;
};
} // namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFFTWComplexToComplexImageFilter.hxx"
#endif
#endif //__itkFFTWComplexToComplexImageFilter_h
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkFFTWComplexToComplexImageFilter_hxx
#define __itkFFTWComplexToComplexImageFilter_hxx
#include "itkFFTWComplexToComplexImageFilter.h"
#include <iostream>
#include "itkIndent.h"
#include "itkMetaDataObject.h"
#include "itkImageRegionIterator.h"
#include "itkProgressReporter.h"
/*
*
* This code was contributed in the Insight Journal paper:
* "FFT Complex to Complex filters and helper classes"
* by Warfield S.
* http://hdl.handle.net/1926/326
* http://www.insight-journal.org/browse/publication/128
*
*/
namespace itk
{
template< class TImage >
void
FFTWComplexToComplexImageFilter< TImage >::
BeforeThreadedGenerateData()
{
// get pointers to the input and output
typename InputImageType::ConstPointer inputPtr = this->GetInput();
typename OutputImageType::Pointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// we don't have a nice progress to report, but at least this simple line
// reports the beginning and the end of the process
ProgressReporter progress(this, 0, 1);
// allocate output buffer memory
outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
outputPtr->Allocate();
const typename InputImageType::SizeType & outputSize =
outputPtr->GetLargestPossibleRegion().GetSize();
const typename OutputImageType::SizeType & inputSize =
inputPtr->GetLargestPossibleRegion().GetSize();
// figure out sizes
// size of input and output aren't the same which is handled in the superclass,
// sort of.
// the input size and output size only differ in the fastest moving dimension
unsigned int total_outputSize = 1;
unsigned int total_inputSize = 1;
for ( unsigned i = 0; i < ImageDimension; i++ )
{
total_outputSize *= outputSize[i];
total_inputSize *= inputSize[i];
}
int transformDirection = 1;
if ( this->GetTransformDirection() == Superclass::INVERSE )
{
transformDirection = -1;
}
typename FFTWProxyType::PlanType plan;
typename FFTWProxyType::ComplexType * in = (typename FFTWProxyType::ComplexType*) inputPtr->GetBufferPointer();
typename FFTWProxyType::ComplexType * out = (typename FFTWProxyType::ComplexType*) outputPtr->GetBufferPointer();
int flags = m_PlanRigor;
if( !m_CanUseDestructiveAlgorithm )
{
// if the input is about to be destroyed, there is no need to force fftw
// to use an non destructive algorithm. If it is not released however,
// we must be careful to not destroy it.
flags = flags | FFTW_PRESERVE_INPUT;
}
int *sizes = new int[ImageDimension];
for(unsigned int i = 0; i < ImageDimension; i++)
{
sizes[(ImageDimension - 1) - i] = inputSize[i];
}
plan = FFTWProxyType::Plan_dft(ImageDimension,sizes,
in,
out,
transformDirection,
flags,
this->GetNumberOfThreads());
delete [] sizes;
FFTWProxyType::Execute(plan);
FFTWProxyType::DestroyPlan(plan);
}
template <class TImage>
void
FFTWComplexToComplexImageFilter< TImage >::
ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType itkNotUsed(threadId) )
{
//
// Normalize the output if backward transform
//
if ( this->GetTransformDirection() == Superclass::INVERSE )
{
typedef ImageRegionIterator< OutputImageType > IteratorType;
unsigned long total_outputSize = this->GetOutput()->GetRequestedRegion().GetNumberOfPixels();
IteratorType it(this->GetOutput(), outputRegionForThread);
while( !it.IsAtEnd() )
{
PixelType val = it.Value();
val /= total_outputSize;
it.Set(val);
++it;
}
}
}
template< class TImage >
bool
FFTWComplexToComplexImageFilter< TImage >::FullMatrix()
{
return false;
}
template< class TImage >
void
FFTWComplexToComplexImageFilter< TImage >::
UpdateOutputData(DataObject * output)
{
// we need to catch that information now, because it is changed later
// during the pipeline execution, and thus can't be grabbed in
// GenerateData().
m_CanUseDestructiveAlgorithm = this->GetInput()->GetReleaseDataFlag();
Superclass::UpdateOutputData( output );
}
template< class TImage >
void
FFTWComplexToComplexImageFilter< TImage >
::PrintSelf(std::ostream & os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "PlanRigor: " << FFTWGlobalConfiguration::GetPlanRigorName(m_PlanRigor) << " (" << m_PlanRigor << ")" << std::endl;
}
} // namespace itk
#endif // _itkFFTWComplexToComplexImageFilter_hxx
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt 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.
=========================================================================*/
#include "itkMacro.h"
#include "otbImage.h"
#include "itkFFTComplexToComplexImageFilter.h"
int main(int argc, char * argv[])
{
typedef otb::Image< std::complex<float>, 2 > ImageType;
typedef itk::FFTComplexToComplexImageFilter <ImageType> FFTFilterType;
FFTFilterType::Pointer fftFilter = FFTFilterType::New();
fftFilter->DebugOn();
return EXIT_SUCCESS;
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment