size_t localsize[2];
String kernel_name;
+ bool is1d = (flags & DFT_ROWS) != 0 || dft_size == 1;
+ String options = buildOptions;
+
if (rows)
{
globalsize[0] = thread_count; globalsize[1] = dft_size;
localsize[0] = thread_count; localsize[1] = 1;
kernel_name = "fft_multi_radix_rows";
+ if (is1d && (flags & DFT_SCALE))
+ options += " -D DFT_SCALE";
}
else
{
globalsize[0] = dft_size; globalsize[1] = thread_count;
localsize[0] = 1; localsize[1] = thread_count;
kernel_name = "fft_multi_radix_cols";
+ if (flags & DFT_SCALE)
+ options += " -D DFT_SCALE";
}
-
- bool is1d = (flags & DFT_ROWS) != 0 || dft_size == 1;
- String options = buildOptions;
+
if (src.channels() == 1)
options += " -D REAL_INPUT";
if (dst.channels() == 1)
options += " -D CCS_OUTPUT";
if ((is1d && src.channels() == 1) || (rows && (flags & DFT_REAL_OUTPUT)))
options += " -D NO_CONJUGATE";
+ if (is1d)
+ options += " -D IS_1D";
ocl::Kernel k(kernel_name.c_str(), ocl::core::fft_oclsrc, options);
if (k.empty())
barrier(CLK_LOCAL_MEM_FENCE);
}
+#ifdef DFT_SCALE
+#define VAL(x, scale) x*scale
+#else
+#define VAL(x, scale) x
+#endif
+
__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,
__constant float2 * twiddles_ptr, const int t, const int nz)
__constant const float2* twiddles = (__constant float2*) twiddles_ptr;
const int ind = x;
const int block_size = LOCAL_SIZE/kercn;
+#ifdef IS_1D
+ float scale = 1.f/dst_cols;
+#else
+ float scale = 1.f/(dst_cols*dst_rows);
+#endif
#ifndef REAL_INPUT
__global const float2* src = (__global const float2*)(src_ptr + mad24(y, src_step, mad24(x, (int)(sizeof(float)*2), src_offset)));
__global float2* dst = (__global float2*)(dst_ptr + mad24(y, dst_step, dst_offset));
#pragma unroll
for (int i=x; i<cols; i+=block_size)
- dst[i] = smem[i];
+ dst[i] = VAL(smem[i], scale);
#else
// pack row to CCS
__local float* smem_1cn = (__local float*) smem;
__global float* dst = (__global float*)(dst_ptr + mad24(y, dst_step, dst_offset));
for (int i=x; i<dst_cols-1; i+=block_size)
- dst[i+1] = smem_1cn[i+2];
+ dst[i+1] = VAL(smem_1cn[i+2], scale);
if (x == 0)
- dst[0] = smem_1cn[0];
+ dst[0] = VAL(smem_1cn[0], scale);
#endif
}
}
__constant const float2* twiddles = (__constant float2*) twiddles_ptr;
const int ind = y;
const int block_size = LOCAL_SIZE/kercn;
+ float scale = 1.f/(dst_rows*dst_cols);
+
#pragma unroll
for (int i=0; i<kercn; i++)
smem[y+i*block_size] = *((__global const float2*)(src + i*block_size*src_step));
__global uchar* dst = dst_ptr + mad24(y, dst_step, mad24(x, (int)(sizeof(float)*2), dst_offset));
#pragma unroll
for (int i=0; i<kercn; i++)
- *((__global float2*)(dst + i*block_size*dst_step)) = smem[y + i*block_size];
+ *((__global float2*)(dst + i*block_size*dst_step)) = VAL(smem[y + i*block_size], scale);
#else
if (x == 0)
{
__local float* smem_1cn = (__local float*) smem;
__global uchar* dst = dst_ptr + mad24(y+1, dst_step, dst_offset);
for (int i=y; i<dst_rows-1; i+=block_size, dst+=dst_step*block_size)
- *((__global float*) dst) = smem_1cn[i+2];
+ *((__global float*) dst) = VAL(smem_1cn[i+2], scale);
if (y == 0)
- *((__global float*) (dst_ptr + dst_offset)) = smem_1cn[0];
+ *((__global float*) (dst_ptr + dst_offset)) = VAL(smem_1cn[0], scale);
}
else if (x == (dst_cols+1)/2)
{
__local float* smem_1cn = (__local float*) smem;
__global uchar* dst = dst_ptr + mad24(dst_cols-1, (int)sizeof(float), mad24(y+1, dst_step, dst_offset));
for (int i=y; i<dst_rows-1; i+=block_size, dst+=dst_step*block_size)
- *((__global float*) dst) = smem_1cn[i+2];
+ *((__global float*) dst) = VAL(smem_1cn[i+2], scale);
if (y == 0)
- *((__global float*) (dst_ptr + mad24(dst_cols-1, (int)sizeof(float), dst_offset))) = smem_1cn[0];
+ *((__global float*) (dst_ptr + mad24(dst_cols-1, (int)sizeof(float), dst_offset))) = VAL(smem_1cn[0], scale);
}
else
{
__global uchar* dst = dst_ptr + mad24(x, (int)sizeof(float)*2, mad24(y, dst_step, dst_offset - (int)sizeof(float)));
#pragma unroll
for (int i=y; i<dst_rows; i+=block_size, dst+=block_size*dst_step)
- vstore2(smem[i], 0, (__global float*) dst);
+ vstore2(VAL(smem[i], scale), 0, (__global float*) dst);
}
#endif
}
////////////////////////////////////////////////////////////////////////////
// Dft
-PARAM_TEST_CASE(Dft, cv::Size, OCL_FFT_TYPE, bool, bool)
+PARAM_TEST_CASE(Dft, cv::Size, OCL_FFT_TYPE, bool, bool, bool)
{
cv::Size dft_size;
int dft_flags, depth, cn, dft_type;
}
if (GET_PARAM(2))
- dft_flags |= cv::DFT_ROWS;
- //if (GET_PARAM(3))
- // if (dft_type == C2C) dft_flags |= cv::DFT_INVERSE;
- //if (GET_PARAM(3))
- // dft_flags |= cv::DFT_SCALE;
-
- inplace = GET_PARAM(3);
- if (inplace && dft_type == 0)
- inplace = 0;
+ dft_flags |= cv::DFT_ROWS;
+ if (GET_PARAM(3))
+ dft_flags |= cv::DFT_SCALE;
+ //if (GET_PARAM(4))
+ // dft_flags |= cv::DFT_INVERSE;
+ inplace = GET_PARAM(4);
+
+
is1d = (dft_flags & DFT_ROWS) != 0 || dft_size.height == 1;
}
udst = udst(cv::Range(0, udst.rows), cv::Range(0, udst.cols/2 + 1));
}
- Mat gpu = udst.getMat(ACCESS_READ);
+ //Mat gpu = udst.getMat(ACCESS_READ);
//std::cout << src << std::endl;
//std::cout << dst << std::endl;
//std::cout << gpu << std::endl;
cv::Size(512, 1), cv::Size(1280, 768)),
Values((OCL_FFT_TYPE) R2C, (OCL_FFT_TYPE) C2C, (OCL_FFT_TYPE) R2R, (OCL_FFT_TYPE) C2R),
Bool(), // DFT_ROWS
+ Bool(), // DFT_SCALE
Bool() // inplace
)
);
#include "opencv2/highgui.hpp"
#include <stdio.h>
-#include <iostream>
-#include <chrono>
using namespace cv;
using namespace std;
int main(int argc, const char ** argv)
{
- //int cols = 4;
- //int rows = 768;
- //srand(0);
- //Mat input(Size(cols, rows), CV_32FC2);
- //for (int i=0; i<cols; i++)
- // for (int j=0; j<rows; j++)
- // input.at<Vec2f>(j,i) = Vec2f((float) rand() / RAND_MAX, (float) rand() / RAND_MAX);
- //Mat dst;
- //
- //UMat gpu_input, gpu_dst;
- //input.copyTo(gpu_input);
- //auto start = std::chrono::system_clock::now();
- //dft(input, dst, DFT_ROWS);
- //auto cpu_duration = chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start);
- //
- //start = std::chrono::system_clock::now();
- //dft(gpu_input, gpu_dst, DFT_ROWS);
- //auto gpu_duration = chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start);
-
- //double n = norm(dst, gpu_dst);
- //cout << "norm = " << n << endl;
- //cout << "CPU time: " << cpu_duration.count() << "ms" << endl;
- //cout << "GPU time: " << gpu_duration.count() << "ms" << endl;
-
-
help();
CommandLineParser parser(argc, argv, keys);
string filename = parser.get<string>(0);
printf("Cannot read image file: %s\n", filename.c_str());
return -1;
}
-
- Mat small_img = img(Rect(0,0,6,6));
-
- int M = getOptimalDFTSize( small_img.rows );
- int N = getOptimalDFTSize( small_img.cols );
+ int M = getOptimalDFTSize( img.rows );
+ int N = getOptimalDFTSize( img.cols );
Mat padded;
- copyMakeBorder(small_img, padded, 0, M - small_img.rows, 0, N - small_img.cols, BORDER_CONSTANT, Scalar::all(0));
+ copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));
- Mat planes[] = {Mat_<float>(padded), Mat::ones(padded.size(), CV_32F)};
- Mat complexImg, complexImg1, complexInput;
+ Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
+ Mat complexImg;
merge(planes, 2, complexImg);
- Mat realInput;
- padded.convertTo(realInput, CV_32F);
- complexInput = complexImg;
- //cout << complexImg << endl;
- //dft(complexImg, complexImg, DFT_REAL_OUTPUT);
- //cout << "Complex to Complex" << endl;
- //cout << complexImg << endl;
- cout << "Complex input" << endl << complexInput << endl;
- cout << "Real input" << endl << realInput << endl;
-
- dft(complexInput, complexImg1, DFT_COMPLEX_OUTPUT);
- cout << "Complex to Complex image: " << endl;
- cout << endl << complexImg1 << endl;
-
- Mat realImg1;
- dft(complexInput, realImg1, DFT_REAL_OUTPUT);
- cout << "Complex to Real image: " << endl;
- cout << endl << realImg1 << endl;
-
- Mat realOut;
- dft(complexImg1, realOut, DFT_INVERSE | DFT_COMPLEX_OUTPUT);
- cout << "Complex to Complex (inverse):" << endl;
- cout << realOut << endl;
-
- Mat complexOut;
- dft(realImg1, complexOut, DFT_INVERSE | DFT_REAL_OUTPUT | DFT_SCALE);
- cout << "Complex to Real (inverse):" << endl;
- cout << complexOut << endl;
+ dft(complexImg, complexImg);
// compute log(1 + sqrt(Re(DFT(img))**2 + Im(DFT(img))**2))
split(complexImg, planes);