diff --git a/include/rtkFFTHilbertImageFilter.h b/include/rtkFFTHilbertImageFilter.h new file mode 100644 index 000000000..a64920758 --- /dev/null +++ b/include/rtkFFTHilbertImageFilter.h @@ -0,0 +1,112 @@ +/*========================================================================= + * + * Copyright RTK 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 rtkFFTHilbertImageFilter_h +#define rtkFFTHilbertImageFilter_h + +#include +#include "rtkConfiguration.h" +#include "rtkFFTProjectionsConvolutionImageFilter.h" +#include "rtkMacro.h" + +namespace rtk +{ + +/** \class FFTHilbertImageFilter + * \brief Implements the Hilbert transform. + * + * \test rtkhilbertfiltertest.cxx + * + * \author Aurélien Coussat + * + * \ingroup RTK ImageToImageFilter + */ + +template +class ITK_EXPORT FFTHilbertImageFilter + : public rtk::FFTProjectionsConvolutionImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_ASSIGN(FFTHilbertImageFilter); + + /** Standard class type alias. */ + using Self = FFTHilbertImageFilter; + using Superclass = rtk::FFTProjectionsConvolutionImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Some convenient type alias. */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using FFTPrecisionType = TFFTPrecision; + using IndexType = typename InputImageType::IndexType; + using SizeType = typename InputImageType::SizeType; + + using FFTInputImageType = typename Superclass::FFTInputImageType; + using FFTInputImagePointer = typename FFTInputImageType::Pointer; + using FFTOutputImageType = typename Superclass::FFTOutputImageType; + using FFTOutputImagePointer = typename FFTOutputImageType::Pointer; + + /** Standard New method. */ + itkNewMacro(Self); + + itkGetMacro(PixelShift, double); + // The Set macro is redefined to clear the current FFT kernel when a parameter + // is modified. + virtual void + SetPixelShift(const double _arg) + { + itkDebugMacro("setting PixelShift to " << _arg); + CLANG_PRAGMA_PUSH + CLANG_SUPPRESS_Wfloat_equal if (this->m_PixelShift != _arg) + { + this->m_PixelShift = _arg; + this->Modified(); + this->m_KernelFFT = nullptr; + } + CLANG_PRAGMA_POP + } + + /** Runtime information support. */ + itkTypeMacro(FFTHilbertImageFilter, FFTProjectionsConvolutionImageFilter); + +protected: + FFTHilbertImageFilter() = default; + ~FFTHilbertImageFilter() override = default; + + void + GenerateOutputInformation() override; + + /** Create and return a pointer to one line of the Hilbert kernel in Fourier space. + * Used in generate data functions. */ + void + UpdateFFTProjectionsConvolutionKernel(const SizeType s) override; + +private: + SizeType m_PreviousKernelUpdateSize; + double m_PixelShift; + +}; // end of class + +} // end namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkFFTHilbertImageFilter.hxx" +#endif + +#endif diff --git a/include/rtkFFTHilbertImageFilter.hxx b/include/rtkFFTHilbertImageFilter.hxx new file mode 100644 index 000000000..926b8ec5d --- /dev/null +++ b/include/rtkFFTHilbertImageFilter.hxx @@ -0,0 +1,101 @@ +/*========================================================================= + * + * Copyright RTK 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 rtkFFTHilbertImageFilter_hxx +#define rtkFFTHilbertImageFilter_hxx + +#include "rtkFFTHilbertImageFilter.h" +#include "itkForwardFFTImageFilter.h" + +namespace rtk +{ + +template +void +FFTHilbertImageFilter::GenerateOutputInformation() +{ + FFTProjectionsConvolutionImageFilter::GenerateOutputInformation(); + + auto origin = this->GetOutput()->GetOrigin(); + double spacing = this->GetOutput()->GetSpacing()[0]; + origin[0] += m_PixelShift * spacing; + this->GetOutput()->SetOrigin(origin); +} + +template +void +FFTHilbertImageFilter::UpdateFFTProjectionsConvolutionKernel(const SizeType s) +{ + if (this->m_KernelFFT.GetPointer() != nullptr && s == this->m_PreviousKernelUpdateSize) + { + return; + } + m_PreviousKernelUpdateSize = s; + + const int width = s[0]; + const int height = s[1]; + + // Allocate kernel + SizeType size; + size.Fill(1); + size[0] = width; + FFTInputImagePointer kernel = FFTInputImageType::New(); + kernel->SetRegions(size); + kernel->Allocate(); + kernel->FillBuffer(0.); + + itk::ImageRegionIterator it(kernel, kernel->GetLargestPossibleRegion()); + + // Compute band-limited kernel in space domain + double spacing = this->GetInput()->GetSpacing()[0]; + IndexType ix; + + while (!it.IsAtEnd()) + { + ix = it.GetIndex(); + + double x = (ix[0] + m_PixelShift) * spacing; + if (x > (size[0] / 2) * spacing) + x -= size[0] * spacing; + + if (x == 0.) + { + it.Set(0.); + } + else + { + double v = spacing * (1. - std::cos(itk::Math::pi * x / spacing)) / (itk::Math::pi * x); + it.Set(v); + } + + ++it; + } + + // FFT kernel + using FFTType = itk::RealToHalfHermitianForwardFFTImageFilter; + typename FFTType::Pointer fftK = FFTType::New(); + fftK->SetInput(kernel); + fftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK->Update(); + + this->m_KernelFFT = fftK->GetOutput(); + this->m_KernelFFT->DisconnectPipeline(); +} + +} // end namespace rtk +#endif diff --git a/wrapping/rtkFFTHilbertImageFilter.wrap b/wrapping/rtkFFTHilbertImageFilter.wrap new file mode 100644 index 000000000..9b51c14ab --- /dev/null +++ b/wrapping/rtkFFTHilbertImageFilter.wrap @@ -0,0 +1,12 @@ +if(RTK_USE_CUDA) + itk_wrap_include(itkCudaImage.h) +endif() + +itk_wrap_class("rtk::FFTHilbertImageFilter" POINTER) + foreach(t ${WRAP_ITK_REAL}) + itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3${ITKM_D}" "itk::Image<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 3>, ${ITKT_D}") + endforeach() + if(RTK_USE_CUDA) + itk_wrap_template("CI${ITKM_F}3CI${ITKM_F}3${ITKM_F}" "itk::CudaImage<${ITKT_F}, 3>, itk::CudaImage<${ITKT_F}, 3>, ${ITKT_F}") + endif() +itk_end_wrap_class()