Introduction
Overview
OpenCV, originally developed by Intel, is a widely used cross-platform computer vision library in the world. It was designed fordevel computational efficiency and with a strong focus on real-time applications. OpenCV is written in C++ and its primary interface is in C++, but it still retains a less comprehensive though extensive C interface.
Intel® Math Kernel Library (Intel® MKL) is a library that accelerates mathematical processing and neural network routines to increase performance while reduce development time. It involves Linear Algebra library, Fast Fourier Transforms (FFT), Vector Math and Statistics functions, and so on.
In this article, how to perform FFT with Intel® MKL and OpenCV on an image will be introduced.
Preliminary Requirements
Intel MKL : Please make sure Intel MKL has been installed into your system. Header files and libraries of Intel MKL can be found in /opt/intel/system_studio_<maversion>/mkl, shown as Fig.1. If Intel MKL is installed with Intel® Parallel Studio, they can be found in /opt/intel/compilers_and_libraries_<version>.
Users can download OpenCV source file package from https://opencv.org/ and follow the documentation on https://docs.opencv.org/master/df/d65/tutorial_table_of_content_introduction.html to build the library for their target platforms. After installing OpenCV into its default directory, header files can found in /usr/local/
include/opencv and /usr/local/include/opencv2, while libraries can be found in /usr/local/lib.
Program with Intel® Math Kernel Library and OpenCV
Use OpenCV functions to read image files
OpenCV provides a number of powerful functions which make reading and writing image files easy and convenient. Meanwhile, these functions support major image formats, like jpeg, png, etc. As introduced in the Overview, OpenCV provides both C and C++ interfaces. In this article, functions with C interface is used.
/* load original image */
IplImage *img_src = cvLoadImage("<IMAGE_NAME>", CV_LOAD_IMAGE_GRAYSCALE);
if (img_src == 0)
{
cout <<"cannot load file" << endl;
return 0;
}
Create necessary variables for data storage
To support FFT calculation, several variables are needed to be created. Please refer to comment lines in the following code snippet to find explanations to each variable.
/* Width, height of the source image
step is a variable indicating data size of one row in memory
*/
int width = img_src->width;
int height = img_src->height;
int step = img_src->widthStep;
/* Storage for source data.
MKL FFT calculation involves floating-point data, thus the grayscale image data is needed
to be transformed into corresponding floating-point array.
*/
double *x_real = (double*)mkl_malloc(width*height*sizeof(double), 64);
uchar *img_src_data = (uchar*)img_src->imageData;
for (i = 0; i < height; ++i)
{
for (j = 0; j < width; ++j)
x_real[i*width+j] = (double)img_src_data[i*step+j];
}
/* Storage for raw FFT result
Result of FFT will be saved as an array of complex numbers.
MKL provides a data type called MKL_Complex16. According to the Fourier Transform theory,
FFT result is sized N/2+1 for 1D FFT, and (W/2+1)*H for 2D FFT.
MKL also provides malloc function called mkl_malloc to make sure memory size of the vari-
able is 4K (default memory page size) aligned. This ensures a high-speed calculation.
*/
MKL_Complex16 *x_out = (MKL_Complex16*)mkl_malloc((width/2+1)*height*sizeof(MKL_Complex16), 64);
/* Storage for buffering FFT result
*/
double *x_fft = (double*)mkl_malloc(width*height*sizeof(double), 64);
/* Storage for FFT result image
FFT result image should be grayscale and hold the same size as the original one.
*/
IplImage *img_fft = cvCreateImage(cvSize(img_src->width, img_src->height), IPL_DEPTH_8U, 1);
Initialize IntelMath Kernel Library FFT routines environment
There are 5 steps that users must follow to perform MKL FFT calculation, which can be found at https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/fourier-transform-functions.html.
MKL_LONG status = 0; // variable to retrieve FFT function's status
DFTI_DESCRIPTOR_HANDLE hand = 0; // Handle of FFT context
MKL_LONG N[2]; // array that holds data size
N[0] = height;
N[1] = width;
// Define data precision to be DOUBLE, and in REAL format
DftiCreateDescriptor(&hand, DFTI_DOUBLE, DFTI_REAL, 2, N);
// Use DftiSetValue to set value to predefined variable.
// DFTI_NOT_INPLACE means don't overwrite original data with FFT results
DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
// Define FFT result data storage format
DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
MKL_LONG rs[3]; // set input data strides
rs[0] = 0;
rs[1] = width;
rs[2] = 1;
MKL_LONG cs[3]; // set output data strides
cs[0] = 0;
cs[1] = width/2+1;
cs[2] = 1;
DftiSetValue(hand, DFTI_INPUT_STRIDES, rs);
DftiSetValue(hand, DFTI_OUTPUT_STRIDES, cs);
DftiSetValue(hand, DFTI_FORWARD_SCALE, (double)1.0/(width*height));
DftiCommitDescriptor(hand); // initialize or commit changes to the FFT descriptor
Note: DFTI_INPUT_STRIDES and DFTI_OUTPUT_STRIDES defines the layout of input and output data. Please refer to https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/dfti-input-strides-dfti-output-strides.html for detailed information about how to determine these two variables.
Note: If the layouts in the forward-domain and backward-domain differ, set these two variables explicitly and commit the descriptor before calling computation functions.
Perform Intel® Math Kernel Library FFT routines calculation
DftiComputeForward and DftiComputeBackward are core functions for FFT calculation.
DftiComputeForward(hand, x_real, x_out);
Unpack FFT result and form 2D FFT power spectrum
2D FFT power spectrum shows power spectrum (amplitude of complex number) as a square, center area corresponds to parts with low frequencies, and outer area corresponds to those with high frequencies. 2D FFT power spectrum can not only show frequency features of a 2-dimentional signal, but also bring convenience to image/signal processing, like low-pass or Gaussian filtering/image blur.
In this section, how to unpack the IntelMKL FFT result to form 2D FFT power spectrum is introduced. In this example, packing format of the result is CCE, because DFTI_CONJUGATE_EVEN_STORAGE is set to be DFTI_COMPLEX_COMP
LEX.
double min = 10000000000000000000.0;
double max = 0;
for (i = 0; i < height; i++)
{
for (j = 0; j < width; j++)
{
MKL_Complex16 val;
if(j < width/2+1)
{
val.real = x_out[i*(width/2+1)+j].real;
val.imag = x_out[i*(width/2+1)+j].imag;
}
else // unpack CCE format result
{
if(i == 0)
{
val.real = x_out[width-j].real;
val.imag = x_out[width-j].imag *(-1);
}
else
{
val.real = x_out[(height-i)*(width/2+1)+width-j].real;
val.imag = x_out[(height-i)*(width/2+1)+width-j].imag *(-1);
}
}
double amp = log(sqrt(val.real*val.real+val.imag*val.imag));
x_fft[i*width+j] = amp; // save amplitude of individual complex number into buffer
if(amp < min)
min = amp;
if(amp > max)
max = amp;
}
}
for (i = 0; i < height*width; i++)
x_fft[i] = 255.0 * (x_fft[i] - min) / (double)(max-min); // normalize buffer values to [0-255]
int i0, j0;
uchar* img_fft_data = (uchar*)img_fft->imageData;
for (i = 0; i < height; ++i) // shift zero-frequency component to center
{
for (j = 0; j < width; ++j)
{
if (i < height / 2)
i0 = i + height / 2;
else
i0 = i - height / 2;
if (j < width / 2)
j0 = j + width / 2;
else
j0 = j - width / 2;
img_fft_data[i * step + j] = (uchar)x_fft[i0 * width + j0];
}
}
Data stored in img_fft is the 2D FFT power spectrum of the source image. It can be displayed or written into disk using OpenCV functions.
Release resources
Manually allocated memory should be released if there are no further usage.
cvReleaseImage(&img_src);
cvReleaseImage(&img_fft);
DftiFreeDescriptor(&hand);
mkl_free(x_real);
mkl_free(x_fft);
mkl_free(x_out);
Results
Figure [img:results] illustrates result of this example. Image in the left is the source image, while the one in the right is the 2D FFT power spectrum.
A source image and the 2D FFT power spectrum