String buildOptions;
int thread_count;
int dft_size;
+ int dft_depth;
bool status;
public:
- OCL_FftPlan(int _size) : dft_size(_size), status(true)
+ OCL_FftPlan(int _size, int _depth) : dft_size(_size), dft_depth(_depth), status(true)
{
+ CV_Assert( dft_depth == CV_32F || dft_depth == CV_64F );
+
int min_radix;
std::vector<int> radixes, blocks;
ocl_getRadixes(dft_size, radixes, blocks, min_radix);
n *= radix;
}
- twiddles.create(1, twiddle_size, CV_32FC2);
- Mat tw = twiddles.getMat(ACCESS_WRITE);
- float* ptr = tw.ptr<float>();
- int ptr_index = 0;
-
- n = 1;
- for (size_t i=0; i<radixes.size(); i++)
- {
- int radix = radixes[i];
- n *= radix;
-
- for (int j=1; j<radix; j++)
- {
- double theta = -CV_2PI*j/n;
-
- for (int k=0; k<(n/radix); k++)
- {
- ptr[ptr_index++] = (float) cos(k*theta);
- ptr[ptr_index++] = (float) sin(k*theta);
- }
- }
- }
+ twiddles.create(1, twiddle_size, CV_MAKE_TYPE(dft_depth, 2));
+ if (dft_depth == CV_32F)
+ fillRadixTable<float>(twiddles, radixes);
+ else
+ fillRadixTable<double>(twiddles, radixes);
- buildOptions = format("-D LOCAL_SIZE=%d -D kercn=%d -D RADIX_PROCESS=%s",
- dft_size, min_radix, radix_processing.c_str());
+ buildOptions = format("-D LOCAL_SIZE=%d -D kercn=%d -D FT=%s -D CT=%s%s -D RADIX_PROCESS=%s",
+ dft_size, min_radix, ocl::typeToStr(dft_depth), ocl::typeToStr(CV_MAKE_TYPE(dft_depth, 2)),
+ dft_depth == CV_64F ? " -D DOUBLE_SUPPORT" : "", radix_processing.c_str());
}
bool enqueueTransform(InputArray _src, OutputArray _dst, int num_dfts, int flags, int fftType, bool rows = true) const
UMat src = _src.getUMat();
UMat dst = _dst.getUMat();
- int type = src.type(), depth = CV_MAT_DEPTH(type);
-
size_t globalsize[2];
size_t localsize[2];
String kernel_name;
bool is1d = (flags & DFT_ROWS) != 0 || num_dfts == 1;
bool inv = (flags & DFT_INVERSE) != 0;
- String options = buildOptions + format(" -D FT=%s CT=%s%s", ocl::typeToStr(depth),
- ocl::typeToStr(CV_MAKE_TYPE(depth, 2)),
- depth == CV_64F ? " -D DOUBLE_SUPPORT" : "");
+ String options = buildOptions;
if (rows)
{
if (k.empty())
return false;
- k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst), ocl::KernelArg::PtrReadOnly(twiddles), thread_count, num_dfts);
+ k.args(ocl::KernelArg::ReadOnly(src), ocl::KernelArg::WriteOnly(dst), ocl::KernelArg::ReadOnlyNoSize(twiddles), thread_count, num_dfts);
return k.run(2, globalsize, localsize, false);
}
min_radix = min(min_radix, block*radix);
}
}
+
+ template <typename T>
+ static void fillRadixTable(UMat twiddles, const std::vector<int>& radixes)
+ {
+ Mat tw = twiddles.getMat(ACCESS_WRITE);
+ T* ptr = tw.ptr<T>();
+ int ptr_index = 0;
+
+ int n = 1;
+ for (size_t i=0; i<radixes.size(); i++)
+ {
+ int radix = radixes[i];
+ n *= radix;
+
+ for (int j=1; j<radix; j++)
+ {
+ double theta = -CV_2PI*j/n;
+
+ for (int k=0; k<(n/radix); k++)
+ {
+ ptr[ptr_index++] = (T) cos(k*theta);
+ ptr[ptr_index++] = (T) sin(k*theta);
+ }
+ }
+ }
+ }
};
class OCL_FftPlanCache
return planCache;
}
- Ptr<OCL_FftPlan> getFftPlan(int dft_size)
+ Ptr<OCL_FftPlan> getFftPlan(int dft_size, int depth)
{
- std::map<int, Ptr<OCL_FftPlan> >::iterator f = planStorage.find(dft_size);
+ int key = (dft_size << 16) | (depth & 0xFFFF);
+ std::map<int, Ptr<OCL_FftPlan> >::iterator f = planStorage.find(key);
if (f != planStorage.end())
{
return f->second;
}
else
{
- Ptr<OCL_FftPlan> newPlan = Ptr<OCL_FftPlan>(new OCL_FftPlan(dft_size));
- planStorage[dft_size] = newPlan;
+ Ptr<OCL_FftPlan> newPlan = Ptr<OCL_FftPlan>(new OCL_FftPlan(dft_size, depth));
+ planStorage[key] = newPlan;
return newPlan;
}
}
static bool ocl_dft_rows(InputArray _src, OutputArray _dst, int nonzero_rows, int flags, int fftType)
{
- Ptr<OCL_FftPlan> plan = OCL_FftPlanCache::getInstance().getFftPlan(_src.cols());
+ int type = _src.type(), depth = CV_MAT_DEPTH(type);
+ Ptr<OCL_FftPlan> plan = OCL_FftPlanCache::getInstance().getFftPlan(_src.cols(), depth);
return plan->enqueueTransform(_src, _dst, nonzero_rows, flags, fftType, true);
}
static bool ocl_dft_cols(InputArray _src, OutputArray _dst, int nonzero_cols, int flags, int fftType)
{
- Ptr<OCL_FftPlan> plan = OCL_FftPlanCache::getInstance().getFftPlan(_src.rows());
+ int type = _src.type(), depth = CV_MAT_DEPTH(type);
+ Ptr<OCL_FftPlan> plan = OCL_FftPlanCache::getInstance().getFftPlan(_src.rows(), depth);
return plan->enqueueTransform(_src, _dst, nonzero_cols, flags, fftType, false);
}
{
int type = _src.type(), cn = CV_MAT_CN(type), depth = CV_MAT_DEPTH(type);
Size ssize = _src.size();
- bool doubleSupport = ocl::Device::getDefault().doubleFPConfig();
+ bool doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0;
if ( !((cn == 1 || cn == 2) && (depth == CV_32F || (depth == CV_64F && doubleSupport))) )
return false;
#endif
__attribute__((always_inline))
-float2 mul_float2(float2 a, float2 b) {
- return (float2)(fma(a.x, b.x, -a.y * b.y), fma(a.x, b.y, a.y * b.x));
+CT mul_complex(CT a, CT b) {
+ return (CT)(fma(a.x, b.x, -a.y * b.y), fma(a.x, b.y, a.y * b.x));
}
__attribute__((always_inline))
-float2 twiddle(float2 a) {
- return (float2)(a.y, -a.x);
+CT twiddle(CT a) {
+ return (CT)(a.y, -a.x);
}
__attribute__((always_inline))
-void butterfly2(float2 a0, float2 a1, __local float2* smem, __global const float2* twiddles,
+void butterfly2(CT a0, CT a1, __local CT* smem, __global const CT* twiddles,
const int x, const int block_size)
{
const int k = x & (block_size - 1);
- a1 = mul_float2(twiddles[k], a1);
+ a1 = mul_complex(twiddles[k], a1);
const int dst_ind = (x << 1) - k;
smem[dst_ind] = a0 + a1;
}
__attribute__((always_inline))
-void butterfly4(float2 a0, float2 a1, float2 a2, float2 a3, __local float2* smem, __global const float2* twiddles,
+void butterfly4(CT a0, CT a1, CT a2, CT a3, __local CT* smem, __global const CT* twiddles,
const int x, const int block_size)
{
const int k = x & (block_size - 1);
- a1 = mul_float2(twiddles[k], a1);
- a2 = mul_float2(twiddles[k + block_size], a2);
- a3 = mul_float2(twiddles[k + 2*block_size], a3);
+ a1 = mul_complex(twiddles[k], a1);
+ a2 = mul_complex(twiddles[k + block_size], a2);
+ a3 = mul_complex(twiddles[k + 2*block_size], a3);
const int dst_ind = ((x - k) << 2) + k;
- float2 b0 = a0 + a2;
+ CT b0 = a0 + a2;
a2 = a0 - a2;
- float2 b1 = a1 + a3;
+ CT b1 = a1 + a3;
a3 = twiddle(a1 - a3);
smem[dst_ind] = b0 + b1;
}
__attribute__((always_inline))
-void butterfly3(float2 a0, float2 a1, float2 a2, __local float2* smem, __global const float2* twiddles,
+void butterfly3(CT a0, CT a1, CT a2, __local CT* smem, __global const CT* twiddles,
const int x, const int block_size)
{
const int k = x % block_size;
- a1 = mul_float2(twiddles[k], a1);
- a2 = mul_float2(twiddles[k+block_size], a2);
+ a1 = mul_complex(twiddles[k], a1);
+ a2 = mul_complex(twiddles[k+block_size], a2);
const int dst_ind = ((x - k) * 3) + k;
- float2 b1 = a1 + a2;
+ CT b1 = a1 + a2;
a2 = twiddle(sin_120*(a1 - a2));
- float2 b0 = a0 - (float2)(0.5f)*b1;
+ CT b0 = a0 - (CT)(0.5f)*b1;
smem[dst_ind] = a0 + b1;
smem[dst_ind + block_size] = b0 + a2;
}
__attribute__((always_inline))
-void butterfly5(float2 a0, float2 a1, float2 a2, float2 a3, float2 a4, __local float2* smem, __global const float2* twiddles,
+void butterfly5(CT a0, CT a1, CT a2, CT a3, CT a4, __local CT* smem, __global const CT* twiddles,
const int x, const int block_size)
{
const int k = x % block_size;
- a1 = mul_float2(twiddles[k], a1);
- a2 = mul_float2(twiddles[k + block_size], a2);
- a3 = mul_float2(twiddles[k+2*block_size], a3);
- a4 = mul_float2(twiddles[k+3*block_size], a4);
+ a1 = mul_complex(twiddles[k], a1);
+ a2 = mul_complex(twiddles[k + block_size], a2);
+ a3 = mul_complex(twiddles[k+2*block_size], a3);
+ a4 = mul_complex(twiddles[k+3*block_size], a4);
const int dst_ind = ((x - k) * 5) + k;
- __local float2* dst = smem + dst_ind;
+ __local CT* dst = smem + dst_ind;
- float2 b0, b1, b5;
+ CT b0, b1, b5;
b1 = a1 + a4;
a1 -= a4;
a3 -= a2;
a2 = b1 + a4;
- b0 = a0 - (float2)0.25f * a2;
+ b0 = a0 - (CT)0.25f * a2;
b1 = fft5_2 * (b1 - a4);
- a4 = fft5_3 * (float2)(-a1.y - a3.y, a1.x + a3.x);
- b5 = (float2)(a4.x - fft5_5 * a1.y, a4.y + fft5_5 * a1.x);
+ a4 = fft5_3 * (CT)(-a1.y - a3.y, a1.x + a3.x);
+ b5 = (CT)(a4.x - fft5_5 * a1.y, a4.y + fft5_5 * a1.x);
a4.x += fft5_4 * a3.y;
a4.y -= fft5_4 * a3.x;
}
__attribute__((always_inline))
-void fft_radix2(__local float2* smem, __global const float2* twiddles, const int x, const int block_size, const int t)
+void fft_radix2(__local CT* smem, __global const CT* twiddles, const int x, const int block_size, const int t)
{
- float2 a0, a1;
+ CT a0, a1;
if (x < t)
{
}
__attribute__((always_inline))
-void fft_radix2_B2(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix2_B2(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/2;
- float2 a0, a1, a2, a3;
+ CT a0, a1, a2, a3;
if (x1 < t/2)
{
}
__attribute__((always_inline))
-void fft_radix2_B3(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix2_B3(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/3;
const int x3 = x1 + 2*t/3;
- float2 a0, a1, a2, a3, a4, a5;
+ CT a0, a1, a2, a3, a4, a5;
if (x1 < t/3)
{
}
__attribute__((always_inline))
-void fft_radix2_B4(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix2_B4(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int thread_block = t/4;
const int x2 = x1 + thread_block;
const int x3 = x1 + 2*thread_block;
const int x4 = x1 + 3*thread_block;
- float2 a0, a1, a2, a3, a4, a5, a6, a7;
+ CT a0, a1, a2, a3, a4, a5, a6, a7;
if (x1 < t/4)
{
}
__attribute__((always_inline))
-void fft_radix2_B5(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix2_B5(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int thread_block = t/5;
const int x2 = x1 + thread_block;
const int x3 = x1 + 2*thread_block;
const int x4 = x1 + 3*thread_block;
const int x5 = x1 + 4*thread_block;
- float2 a0, a1, a2, a3, a4, a5, a6, a7, a8, a9;
+ CT a0, a1, a2, a3, a4, a5, a6, a7, a8, a9;
if (x1 < t/5)
{
}
__attribute__((always_inline))
-void fft_radix4(__local float2* smem, __global const float2* twiddles, const int x, const int block_size, const int t)
+void fft_radix4(__local CT* smem, __global const CT* twiddles, const int x, const int block_size, const int t)
{
- float2 a0, a1, a2, a3;
+ CT a0, a1, a2, a3;
if (x < t)
{
}
__attribute__((always_inline))
-void fft_radix4_B2(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix4_B2(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/2;
- float2 a0, a1, a2, a3, a4, a5, a6, a7;
+ CT a0, a1, a2, a3, a4, a5, a6, a7;
if (x1 < t/2)
{
}
__attribute__((always_inline))
-void fft_radix4_B3(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix4_B3(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/3;
const int x3 = x2 + t/3;
- float2 a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11;
+ CT a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11;
if (x1 < t/3)
{
}
__attribute__((always_inline))
-void fft_radix8(__local float2* smem, __global const float2* twiddles, const int x, const int block_size, const int t)
+void fft_radix8(__local CT* smem, __global const CT* twiddles, const int x, const int block_size, const int t)
{
const int k = x % block_size;
- float2 a0, a1, a2, a3, a4, a5, a6, a7;
+ CT a0, a1, a2, a3, a4, a5, a6, a7;
if (x < t)
{
int tw_ind = block_size / 8;
a0 = smem[x];
- a1 = mul_float2(twiddles[k], smem[x + t]);
- a2 = mul_float2(twiddles[k + block_size],smem[x+2*t]);
- a3 = mul_float2(twiddles[k+2*block_size],smem[x+3*t]);
- a4 = mul_float2(twiddles[k+3*block_size],smem[x+4*t]);
- a5 = mul_float2(twiddles[k+4*block_size],smem[x+5*t]);
- a6 = mul_float2(twiddles[k+5*block_size],smem[x+6*t]);
- a7 = mul_float2(twiddles[k+6*block_size],smem[x+7*t]);
+ a1 = mul_complex(twiddles[k], smem[x + t]);
+ a2 = mul_complex(twiddles[k + block_size],smem[x+2*t]);
+ a3 = mul_complex(twiddles[k+2*block_size],smem[x+3*t]);
+ a4 = mul_complex(twiddles[k+3*block_size],smem[x+4*t]);
+ a5 = mul_complex(twiddles[k+4*block_size],smem[x+5*t]);
+ a6 = mul_complex(twiddles[k+5*block_size],smem[x+6*t]);
+ a7 = mul_complex(twiddles[k+6*block_size],smem[x+7*t]);
- float2 b0, b1, b6, b7;
+ CT b0, b1, b6, b7;
b0 = a0 + a4;
a4 = a0 - a4;
b1 = a1 + a5;
a5 = a1 - a5;
- a5 = (float2)(SQRT_2) * (float2)(a5.x + a5.y, -a5.x + a5.y);
+ a5 = (CT)(SQRT_2) * (CT)(a5.x + a5.y, -a5.x + a5.y);
b6 = twiddle(a2 - a6);
a2 = a2 + a6;
b7 = a3 - a7;
- b7 = (float2)(SQRT_2) * (float2)(-b7.x + b7.y, -b7.x - b7.y);
+ b7 = (CT)(SQRT_2) * (CT)(-b7.x + b7.y, -b7.x - b7.y);
a3 = a3 + a7;
a0 = b0 + a2;
if (x < t)
{
const int dst_ind = ((x - k) << 3) + k;
- __local float2* dst = smem + dst_ind;
+ __local CT* dst = smem + dst_ind;
dst[0] = a0 + a1;
dst[block_size] = a4 + a5;
}
__attribute__((always_inline))
-void fft_radix3(__local float2* smem, __global const float2* twiddles, const int x, const int block_size, const int t)
+void fft_radix3(__local CT* smem, __global const CT* twiddles, const int x, const int block_size, const int t)
{
- float2 a0, a1, a2;
+ CT a0, a1, a2;
if (x < t)
{
}
__attribute__((always_inline))
-void fft_radix3_B2(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix3_B2(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/2;
- float2 a0, a1, a2, a3, a4, a5;
+ CT a0, a1, a2, a3, a4, a5;
if (x1 < t/2)
{
}
__attribute__((always_inline))
-void fft_radix3_B3(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix3_B3(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1 + t/3;
const int x3 = x2 + t/3;
- float2 a0, a1, a2, a3, a4, a5, a6, a7, a8;
+ CT a0, a1, a2, a3, a4, a5, a6, a7, a8;
if (x1 < t/3)
{
}
__attribute__((always_inline))
-void fft_radix3_B4(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix3_B4(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int thread_block = t/4;
const int x2 = x1 + thread_block;
const int x3 = x1 + 2*thread_block;
const int x4 = x1 + 3*thread_block;
- float2 a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11;
+ CT a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11;
if (x1 < t/4)
{
}
__attribute__((always_inline))
-void fft_radix5(__local float2* smem, __global const float2* twiddles, const int x, const int block_size, const int t)
+void fft_radix5(__local CT* smem, __global const CT* twiddles, const int x, const int block_size, const int t)
{
const int k = x % block_size;
- float2 a0, a1, a2, a3, a4;
+ CT a0, a1, a2, a3, a4;
if (x < t)
{
}
__attribute__((always_inline))
-void fft_radix5_B2(__local float2* smem, __global const float2* twiddles, const int x1, const int block_size, const int t)
+void fft_radix5_B2(__local CT* smem, __global const CT* twiddles, const int x1, const int block_size, const int t)
{
const int x2 = x1+t/2;
- float2 a0, a1, a2, a3, a4, a5, a6, a7, a8, a9;
+ CT a0, a1, a2, a3, a4, a5, a6, a7, a8, a9;
if (x1 < t/2)
{
__kernel void fft_multi_radix_rows(__global const uchar* src_ptr, int src_step, int src_offset, int src_rows, int src_cols,
__global uchar* dst_ptr, int dst_step, int dst_offset, int dst_rows, int dst_cols,
- __global float2* twiddles_ptr, const int t, const int nz)
+ __global CT* twiddles_ptr, int twiddles_step, int twiddles_offset, const int t, const int nz)
{
const int x = get_global_id(0);
const int y = get_group_id(1);
if (y < nz)
{
__local CT smem[LOCAL_SIZE];
- __global const float2* twiddles = (__global float2*) twiddles_ptr;
+ __global const CT* twiddles = (__global const CT*)(twiddles_ptr + twiddles_offset);
const int ind = x;
#ifdef IS_1D
FT scale = (FT) 1/dst_cols;
__kernel void fft_multi_radix_cols(__global const uchar* src_ptr, int src_step, int src_offset, int src_rows, int src_cols,
__global uchar* dst_ptr, int dst_step, int dst_offset, int dst_rows, int dst_cols,
- __global float2* twiddles_ptr, const int t, const int nz)
+ __global CT* twiddles_ptr, int twiddles_step, int twiddles_offset, const int t, const int nz)
{
const int x = get_group_id(0);
const int y = get_global_id(1);
{
__local CT smem[LOCAL_SIZE];
__global const uchar* src = src_ptr + mad24(y, src_step, mad24(x, (int)(sizeof(CT)), src_offset));
- __global const float2* twiddles = (__global float2*) twiddles_ptr;
+ __global const CT* twiddles = (__global const CT*)(twiddles_ptr + twiddles_offset);
const int ind = y;
const int block_size = LOCAL_SIZE/kercn;
FT scale = 1.f/(dst_rows*dst_cols);
__kernel void ifft_multi_radix_rows(__global const uchar* src_ptr, int src_step, int src_offset, int src_rows, int src_cols,
__global uchar* dst_ptr, int dst_step, int dst_offset, int dst_rows, int dst_cols,
- __global float2* twiddles_ptr, const int t, const int nz)
+ __global CT* twiddles_ptr, int twiddles_step, int twiddles_offset, const int t, const int nz)
{
const int x = get_global_id(0);
const int y = get_group_id(1);
if (y < nz)
{
__local CT smem[LOCAL_SIZE];
- __global const float2* twiddles = (__global float2*) twiddles_ptr;
+ __global const CT* twiddles = (__global const CT*)(twiddles_ptr + twiddles_offset);
const int ind = x;
#if defined(COMPLEX_INPUT) && !defined(NO_CONJUGATE)
__kernel void ifft_multi_radix_cols(__global const uchar* src_ptr, int src_step, int src_offset, int src_rows, int src_cols,
__global uchar* dst_ptr, int dst_step, int dst_offset, int dst_rows, int dst_cols,
- __global float2* twiddles_ptr, const int t, const int nz)
+ __global CT* twiddles_ptr, int twiddles_step, int twiddles_offset, const int t, const int nz)
{
const int x = get_group_id(0);
const int y = get_global_id(1);
__local CT smem[LOCAL_SIZE];
__global const uchar* src = src_ptr + mad24(y, src_step, mad24(x, (int)(sizeof(CT)), src_offset));
__global uchar* dst = dst_ptr + mad24(y, dst_step, mad24(x, (int)(sizeof(CT)), dst_offset));
- __global const float2* twiddles = (__global float2*) twiddles_ptr;
+ __global const CT* twiddles = (__global const CT*)(twiddles_ptr + twiddles_offset);
const int ind = y;
const int block_size = LOCAL_SIZE/kercn;
#else
if (x < nz)
{
- __global const CT* twiddles = (__global CT*) twiddles_ptr;
+ __global const CT* twiddles = (__global const CT*)(twiddles_ptr + twiddles_offset);
const int ind = y;
const int block_size = LOCAL_SIZE/kercn;