Skip to content

Commit 2674076

Browse files
tbirdsotbirdso
authored
Merge pull request #17 from tbirdso/1d-fft
ENH: Add accelerated 1D FFT classes
2 parents 083e663 + 32f7c2e commit 2674076

21 files changed

Lines changed: 1432 additions & 4 deletions

‎.github/workflows/test-gpu.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ jobs:
1212
opencl-headers-git-tag: "v2021.04.29"
1313
opencl-version: 120
1414
vkfft-backend: 3
15-
itk-git-tag: "c693dbc830859ea10d573301f075e121a4c6c581"
15+
itk-git-tag: "b1a0a91614dc091f9a17c1e4727095c80e9c05ea"
1616
cmake-build-type: "MinSizeRel"
1717
platform-name: "ubuntu-nvidia-gpu"
1818
os: ubuntu-20.04

‎include/itkVkCommon.h‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,9 +63,13 @@ class VkFFTBackend_EXPORT VkCommon
6363

6464
struct VkParameters
6565
{
66-
uint64_t X = 0; // size of fastest varying dimension
67-
uint64_t Y = 1; // size of second-fastest varying dimension, if any, otherwise 1.
68-
uint64_t Z = 1; // size of third-fastest varying dimension, if any, otherwise 1.
66+
uint64_t X = 0; // size of fastest varying dimension
67+
uint64_t Y = 1; // size of second-fastest varying dimension, if any, otherwise 1.
68+
uint64_t Z = 1; // size of third-fastest varying dimension, if any, otherwise 1.
69+
uint64_t omitDimension[3] = { 0,
70+
0,
71+
0 }; // disable FFT for this dimension (0 - FFT enabled, 1 - FFT disabled). Default 0.
72+
// Doesn't work for R2C dimension 0 for now. Doesn't work with convolutions.
6973
PrecisionEnum P = PrecisionEnum::FLOAT; // type for real numbers
7074
uint64_t B = 1; // Number of batches -- always 1
7175
uint64_t N = 1; // Number of redundant iterations, for benchmarking -- always 1.

‎include/itkVkComplexToComplex1DFFTImageFilter.h‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* http://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
#ifndef itkVkComplexToComplex1DFFTImageFilter_h
20+
#define itkVkComplexToComplex1DFFTImageFilter_h
21+
22+
#include "itkVkCommon.h"
23+
#include "itkComplexToComplex1DFFTImageFilter.h"
24+
25+
namespace itk
26+
{
27+
28+
/**
29+
*\class VkComplexToComplex1DFFTImageFilter
30+
*
31+
* \brief Implements an API to enable the 1D Fourier transform or the inverse
32+
* Fourier transform of images with complex valued voxels to be computed using
33+
* the VkFFT library.
34+
*
35+
* This filter is multithreaded and by default supports input images with sizes which are
36+
* divisible by primes up to 13.
37+
*
38+
* Execution on input images with sizes divisible by primes greater than 17 may succeed
39+
* with a fallback on Bluestein's algorithm per VkFFT with a cost to performance and output precision.
40+
*
41+
* \ingroup FourierTransform
42+
* \ingroup MultiThreaded
43+
* \ingroup ITKFFT
44+
* \ingroup VkFFTBackend
45+
*/
46+
template <typename TImage>
47+
class VkComplexToComplex1DFFTImageFilter : public ComplexToComplex1DFFTImageFilter<TImage>
48+
{
49+
public:
50+
ITK_DISALLOW_COPY_AND_MOVE(VkComplexToComplex1DFFTImageFilter);
51+
52+
using InputImageType = TImage;
53+
using OutputImageType = TImage;
54+
static_assert(std::is_same<typename TImage::PixelType, std::complex<float>>::value ||
55+
std::is_same<typename TImage::PixelType, std::complex<double>>::value,
56+
"Unsupported pixel type");
57+
static_assert(TImage::ImageDimension >= 1 && TImage::ImageDimension <= 3, "Unsupported image dimension");
58+
59+
/** Standard class type aliases. */
60+
using Self = VkComplexToComplex1DFFTImageFilter;
61+
using Superclass = ComplexToComplex1DFFTImageFilter<TImage>;
62+
using Pointer = SmartPointer<Self>;
63+
using ConstPointer = SmartPointer<const Self>;
64+
using InputPixelType = typename InputImageType::PixelType;
65+
using OutputPixelType = typename OutputImageType::PixelType;
66+
using ComplexType = InputPixelType;
67+
using RealType = typename ComplexType::value_type;
68+
using SizeType = typename InputImageType::SizeType;
69+
using SizeValueType = typename InputImageType::SizeValueType;
70+
using OutputImageRegionType = typename OutputImageType::RegionType;
71+
72+
/** Method for creation through the object factory. */
73+
itkNewMacro(Self);
74+
75+
/** Run-time type information (and related methods). */
76+
itkTypeMacro(VkComplexToComplex1DFFTImageFilter, ComplexToComplex1DFFTImageFilter);
77+
78+
static constexpr unsigned int ImageDimension = InputImageType::ImageDimension;
79+
80+
itkGetMacro(DeviceID, uint64_t);
81+
itkSetMacro(DeviceID, uint64_t);
82+
83+
SizeValueType
84+
GetSizeGreatestPrimeFactor() const override;
85+
86+
protected:
87+
VkComplexToComplex1DFFTImageFilter() = default;
88+
~VkComplexToComplex1DFFTImageFilter() override = default;
89+
90+
void
91+
GenerateData() override;
92+
93+
void
94+
PrintSelf(std::ostream & os, Indent indent) const override;
95+
96+
private:
97+
uint64_t m_DeviceID{ 0UL };
98+
99+
VkCommon m_VkCommon{};
100+
};
101+
102+
} // namespace itk
103+
104+
#ifndef ITK_MANUAL_INSTANTIATION
105+
# include "itkVkComplexToComplex1DFFTImageFilter.hxx"
106+
#endif
107+
108+
#endif // itkVkComplexToComplex1DFFTImageFilter_h

‎include/itkVkComplexToComplex1DFFTImageFilter.hxx‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* http://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
#ifndef itkVkComplexToComplex1DFFTImageFilter_hxx
20+
#define itkVkComplexToComplex1DFFTImageFilter_hxx
21+
22+
#include "itkVkComplexToComplex1DFFTImageFilter.h"
23+
#include "vkFFT.h"
24+
#include "itkImageRegionIterator.h"
25+
#include "itkIndent.h"
26+
#include "itkMetaDataObject.h"
27+
#include "itkProgressReporter.h"
28+
29+
namespace itk
30+
{
31+
template <typename TImage>
32+
void
33+
VkComplexToComplex1DFFTImageFilter<TImage>::GenerateData()
34+
{
35+
// get pointers to the input and output
36+
const InputImageType * const input{ this->GetInput() };
37+
OutputImageType * const output{ this->GetOutput() };
38+
39+
if (!input || !output)
40+
{
41+
return;
42+
}
43+
44+
// we don't have a nice progress to report, but at least this simple line
45+
// reports the beginning and the end of the process
46+
const ProgressReporter progress(this, 0, 1);
47+
48+
// allocate output buffer memory
49+
output->SetBufferedRegion(output->GetRequestedRegion());
50+
output->Allocate();
51+
52+
const SizeType & inputSize{ input->GetLargestPossibleRegion().GetSize() };
53+
54+
const InputPixelType * const inputCPUBuffer{ input->GetBufferPointer() };
55+
OutputPixelType * const outputCPUBuffer{ output->GetBufferPointer() };
56+
itkAssertOrThrowMacro(inputCPUBuffer != nullptr, "No CPU input buffer");
57+
itkAssertOrThrowMacro(outputCPUBuffer != nullptr, "No CPU output buffer");
58+
const SizeValueType inBytes{ input->GetLargestPossibleRegion().GetNumberOfPixels() * sizeof(InputPixelType) };
59+
const SizeValueType outBytes{ output->GetLargestPossibleRegion().GetNumberOfPixels() * sizeof(OutputPixelType) };
60+
itkAssertOrThrowMacro(inBytes == outBytes, "CPU input and output buffers are of different sizes.");
61+
62+
// Mostly use defaults for VkCommon::VkGPU
63+
typename VkCommon::VkGPU vkGPU;
64+
vkGPU.device_id = m_DeviceID;
65+
66+
// Describe this filter in VkCommon::VkParameters
67+
typename VkCommon::VkParameters vkParameters;
68+
if (ImageDimension > 0)
69+
vkParameters.X = inputSize[0];
70+
if (ImageDimension > 1)
71+
vkParameters.Y = inputSize[1];
72+
if (ImageDimension > 2)
73+
vkParameters.Z = inputSize[2];
74+
if (std::is_same<RealType, float>::value)
75+
vkParameters.P = VkCommon::PrecisionEnum::FLOAT;
76+
else if (std::is_same<RealType, double>::value)
77+
vkParameters.P = VkCommon::PrecisionEnum::DOUBLE;
78+
else
79+
itkAssertOrThrowMacro(false, "Unsupported type for real numbers.");
80+
vkParameters.fft = VkCommon::FFTEnum::C2C;
81+
vkParameters.PSize = sizeof(RealType);
82+
vkParameters.I = this->GetTransformDirection() == Superclass::TransformDirectionType::INVERSE
83+
? VkCommon::DirectionEnum::INVERSE
84+
: VkCommon::DirectionEnum::FORWARD;
85+
vkParameters.normalized = vkParameters.I == VkCommon::DirectionEnum::INVERSE
86+
? VkCommon::NormalizationEnum::NORMALIZED
87+
: VkCommon::NormalizationEnum::UNNORMALIZED;
88+
for (size_t dim = 0; dim < ImageDimension; ++dim)
89+
{
90+
if (this->GetDirection() != dim)
91+
{
92+
vkParameters.omitDimension[dim] = 1; // omit dimensions other than in the given direction.
93+
}
94+
}
95+
96+
vkParameters.inputCPUBuffer = inputCPUBuffer;
97+
vkParameters.inputBufferBytes = inBytes;
98+
vkParameters.outputCPUBuffer = outputCPUBuffer;
99+
vkParameters.outputBufferBytes = outBytes;
100+
101+
const VkFFTResult resFFT{ m_VkCommon.Run(vkGPU, vkParameters) };
102+
if (resFFT != VKFFT_SUCCESS)
103+
{
104+
std::ostringstream mesg;
105+
mesg << "VkFFT third-party library failed with error code " << resFFT << ".";
106+
itkAssertOrThrowMacro(false, mesg.str());
107+
}
108+
}
109+
110+
template <typename TImage>
111+
void
112+
VkComplexToComplex1DFFTImageFilter<TImage>::PrintSelf(std::ostream & os, Indent indent) const
113+
{
114+
Superclass::PrintSelf(os, indent);
115+
os << indent << "DeviceID: " << m_DeviceID << std::endl;
116+
}
117+
118+
template <typename TImage>
119+
typename VkComplexToComplex1DFFTImageFilter<TImage>::SizeValueType
120+
VkComplexToComplex1DFFTImageFilter<TImage>::GetSizeGreatestPrimeFactor() const
121+
{
122+
return m_VkCommon.GetGreatestPrimeFactor();
123+
}
124+
125+
} // end namespace itk
126+
127+
#endif // _itkVkComplexToComplex1DFFTImageFilter_hxx

‎include/itkVkForward1DFFTImageFilter.h‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* http://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkVkForward1DFFTImageFilter_h
19+
#define itkVkForward1DFFTImageFilter_h
20+
21+
#include "itkForward1DFFTImageFilter.h"
22+
#include "itkVkCommon.h"
23+
24+
namespace itk
25+
{
26+
/**
27+
*\class VkForward1DFFTImageFilter
28+
*
29+
* \brief Vk-based forward 1D Fast Fourier Transform.
30+
*
31+
* This filter computes the forward Fourier transform in one dimension of an image. The
32+
* implementation is based on the VkFFT library.
33+
*
34+
* This filter is multithreaded and by default supports input images with sizes which are
35+
* divisible by primes up to 13.
36+
*
37+
* Execution on input images with sizes divisible by primes greater than 17 may succeed
38+
* with a fallback on Bluestein's algorithm per VkFFT with a cost to performance and output precision.
39+
*
40+
* \ingroup FourierTransform
41+
* \ingroup MultiThreaded
42+
* \ingroup ITKFFT
43+
* \ingroup VkFFTBackend
44+
*
45+
* \sa VkGlobalConfiguration
46+
* \sa Forward1DFFTImageFilter
47+
*/
48+
template <typename TInputImage>
49+
class VkForward1DFFTImageFilter
50+
: public Forward1DFFTImageFilter<TInputImage,
51+
Image<std::complex<typename TInputImage::PixelType>, TInputImage::ImageDimension>>
52+
{
53+
public:
54+
ITK_DISALLOW_COPY_AND_MOVE(VkForward1DFFTImageFilter);
55+
56+
using InputImageType = TInputImage;
57+
using OutputImageType = Image<std::complex<typename TInputImage::PixelType>, TInputImage::ImageDimension>;
58+
static_assert(std::is_same<typename TInputImage::PixelType, float>::value ||
59+
std::is_same<typename TInputImage::PixelType, double>::value,
60+
"Unsupported pixel type");
61+
static_assert(TInputImage::ImageDimension >= 1 && TInputImage::ImageDimension <= 3, "Unsupported image dimension");
62+
63+
/** Standard class type aliases. */
64+
using Self = VkForward1DFFTImageFilter;
65+
using Superclass = Forward1DFFTImageFilter<InputImageType, OutputImageType>;
66+
using Pointer = SmartPointer<Self>;
67+
using ConstPointer = SmartPointer<const Self>;
68+
69+
using InputPixelType = typename InputImageType::PixelType;
70+
using OutputPixelType = typename OutputImageType::PixelType;
71+
using ComplexType = OutputPixelType;
72+
using RealType = typename ComplexType::value_type;
73+
using SizeType = typename InputImageType::SizeType;
74+
using SizeValueType = typename InputImageType::SizeValueType;
75+
using OutputImageRegionType = typename OutputImageType::RegionType;
76+
77+
/** Method for creation through the object factory. */
78+
itkNewMacro(Self);
79+
80+
/** Run-time type information (and related methods). */
81+
itkTypeMacro(VkForward1DFFTImageFilter, Forward1DFFTImageFilter);
82+
83+
static constexpr unsigned int ImageDimension = InputImageType::ImageDimension;
84+
85+
itkGetMacro(DeviceID, uint64_t);
86+
itkSetMacro(DeviceID, uint64_t);
87+
88+
SizeValueType
89+
GetSizeGreatestPrimeFactor() const override;
90+
91+
protected:
92+
VkForward1DFFTImageFilter() = default;
93+
~VkForward1DFFTImageFilter() override = default;
94+
95+
void
96+
GenerateData() override;
97+
98+
void
99+
PrintSelf(std::ostream & os, Indent indent) const override;
100+
101+
private:
102+
uint64_t m_DeviceID{ 0UL };
103+
104+
VkCommon m_VkCommon{};
105+
};
106+
107+
} // namespace itk
108+
109+
#ifndef ITK_MANUAL_INSTANTIATION
110+
# include "itkVkForward1DFFTImageFilter.hxx"
111+
#endif
112+
113+
#endif // itkVkForward1DFFTImageFilter_h

0 commit comments

Comments
 (0)