- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
I had a quick look at your code. There are two things you may want to check:
There isn't a s_trig_transform() function defined in the current MKL release. It sould be either s_forward_trig_transform() or s_backward_trig_transform(). Which version of MKL are you using?
For staggered cosine transform, the input array must be size (n+1), where n is the size of transform. And the last element must be set to 0.0. Did you take care of this?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Oh, I accidentatly deleted some lines when sanitizing the code I posted. I'm using a function pointer called s_tring_transform to select forward or backward based on the direction input variable. Here's the lines:
void (*trig_transform)(Float *, DFTI_DESCRIPTOR_HANDLE *, int *, Float *, int *);
if(direction < 0)
trig_transform = s_forward_trig_transform;
else
trig_transform = s_backward_trig_transform;
I also made sure the size of the input array was n+1 and ended with a zero. I've also tried using MKL_COSINE_TRANSFORM and MKL_STAGGERED_COSINE_TRANSFORM, but I still couldn't get the original 2D array back.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I attached a complete example and below the output I get. I print out the original matrix, after the forward_dct, and after the backward_dct. I don't completely understand dct, but I thought I should end up with a matrix similar to what I started with.
11 12 13 14 15 16 17 18
21 22 23 24 25 26 27 28
31 32 33 34 35 36 37 38
41 42 43 44 45 46 47 48
51 52 53 54 55 56 57 58
61 62 63 64 65 66 67 68
71 72 73 74 75 76 77 78
81 82 83 84 85 86 87 88
183.047 4.26973 -13.6921 12.576 -13.3736 13.0362 -13.3722 13.1317
-39.3769 -0.193602 3.04743 -2.74132 2.72896 -2.75311 2.72749 -2.74693
-22.2031 -3.91173 1.05791 -1.36401 1.37637 -1.35223 1.37785 -1.3584
14.6098 3.40551 -0.551685 0.857791 -0.87015 0.846003 -0.871624 0.852176
-22.2031 -3.91173 1.05791 -1.36401 1.37637 -1.35223 1.37785 -1.3584
18.813 3.68572 -0.831899 1.138 -1.15036 1.12622 -1.15184 1.13239
-22.2031 -3.91173 1.05791 -1.36401 1.37637 -1.35223 1.37785 -1.3584
19.7666 3.7493 -0.895476 1.20158 -1.21394 1.18979 -1.21541 1.19597
-6.58086 36.9207 -18.6316 35.9385 -10.7853 49.3869 4.92431 64.7248
20.2915 15.3686 36.3422 20.3508 32.4959 10.9024 20.7864 -0.435492
44.1692 26.1707 32.1184 25.1885 39.9647 38.6369 55.6743 53.9748
34.3433 41.3168 50.394 46.299 46.5477 36.8507 34.8381 25.5127
65.4192 44.9207 53.3684 43.9385 61.2147 57.3869 76.9243 72.7248
55.0329 60.6273 71.0836 65.6095 67.2373 56.1611 55.5277 44.8232
86.6691 63.6707 74.6184 62.6885 82.4647 76.1369 98.1743 91.4748
76.1557 79.5044 92.2064 84.4866 88.3601 75.0382 76.6505 63.7003
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
David,
MKL engineers are going to take a look at the problem you reported. Meanwhile, I'd suggest you to take a look at the linear transform functions in Intel Integrated Performance Primitives (Intel IPP). If what you want is 8x8 DCT, then IPP has already had a function for this. Please take a look at the IPP reference manual, chapter 10: https://software.intel.com/sites/products/documentation/doclib/ipp_sa/71/ipp_manual/ippi.pdf
Thanks.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
[ DCT512x512 image ]
[ Output 512x512 image ]
However, I don't know if the same algorithm is used in MKL but the algorithm in the application I used ( implemented in C/C++ with processing on a GPU / very fast with more than 128 frames ( 512x512 ) per second ) is based on a classic one described on en.wikipedia.org/wiki/Discrete_cosine_transform.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
>>>implemented in C/C++ with processing on a GPU / very fast with more than 128 frames ( 512x512 ) per second ) is based on a classic one described on>>>
What GPU do you have?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Sergey, thanks for the test results. I'm using dct in an advanced image processing algorithm that does a 2D dct on the entire image. I was getting odd results, which is why I tried testing dct on a small 8x8 matrix. At least now I know where the problem is.
It looks like the Intel IPP library would do the 2D transform for me, but I don't have a license for that, just MKL. I wonder how hard this is to implement.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I ended up writing my own 1D dct solution, which was fairly easy, but it is slow. I can at least continue algorithm development now.
Hopefully the MKL will get fixed or I may have to find another library or get an IPP license.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I read this thread with interest some months later and just thought I'd post a folllow-up for anyone that might consider using this library.
For me the attached code reproduces the initial matrix perfectly upon transforming and inverting, using the staggered cosine.
For the staggered sine a warning is issued but the matrix is still reproduced.
I'm using MKL 11.0.5.1
Regards
Rodney
P.S. File attachment doesn't seem to work so here's the code:
#include "cmath"
#include "mkl_trig_transforms.h"
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "mkl.h"
#define _USE_MATH_DEFINES
using namespace std;
void dct2(float** data, MKL_INT num_rows, MKL_INT num_cols, int direction) {
// compute 2D forward cosine transform
DFTI_DESCRIPTOR_HANDLE handle;
MKL_INT stat = 0;
MKL_INT ipar[128];
float * par = new float[(3*num_rows/2) + 2]();
// use function pointers to handle whether float is a float or double
// and direction of the transform
void (*s_trig_transform)(float *, DFTI_DESCRIPTOR_HANDLE *, MKL_INT *, float *, MKL_INT *);
if(direction < 0)
s_trig_transform = s_forward_trig_transform;
else
s_trig_transform = s_backward_trig_transform;
MKL_INT tt_type = MKL_STAGGERED_SINE_TRANSFORM;
// first do the rows
s_init_trig_transform(&num_cols, &tt_type, ipar, par, &stat);
if(stat != 0){
cerr << "Error initializing row transform: " << stat << endl;
return;
}
s_commit_trig_transform(data[0], &handle, ipar, par, &stat);
if(stat != 0){
cerr << "Error committing row transform: " << stat << endl;
return;
}
for(unsigned int row = 0; row < num_rows; ++row){
s_trig_transform(data[row], &handle, ipar, par, &stat);
if(stat != 0){
cerr << "Error in row transform: " << stat << endl;
return;
}
}
free_trig_transform(&handle, ipar, &stat);
if(stat != 0){
cerr << "Error freeing row transform: " << stat << endl;
return;
}
// now do the cols
s_init_trig_transform(&num_rows, &tt_type, ipar, par, &stat);
if(stat != 0){
cerr << "Error initializing col transform: " << stat << endl;
return;
}
float * temp_col = new float[num_rows+1]();
s_commit_trig_transform(temp_col, &handle, ipar, par, &stat);
if(stat != 0){
cerr << "Error committing col transform: " << stat << endl;
return;
}
for(unsigned int col = 0; col < num_cols; ++col){
// copy the column into a contiguous vector
for(unsigned int row = 0; row < num_rows; ++row){
temp_col[row] = data[row][col];
}
s_trig_transform(temp_col, &handle, ipar, par, &stat);
if(stat != 0){
cerr << "Error in col transform: " << stat << endl;
return;
}
// copy the results back
for(unsigned int row = 0; row < num_rows; ++row){
data[row][col] = temp_col[row];
}
}
free_trig_transform(&handle, ipar, &stat);
if(stat != 0){
cerr << "Error freeing col transform: " << stat << endl;
return;
}
delete[] temp_col;
delete[] par;
//#]
}
float log2(float x) {
return log(x)/log(2);
}
int main(int argc, char ** argv){
int size = 8;
float ** temp = new float*[size+1]();
temp[0] = new float[(size+1)*(size+1)]();
for (unsigned int row = 1; row < size+1; row++) {
temp[row] = &(temp[0][row*(size+1)]);
}
for(int i = 0; i < size; ++i){
for(int j = 0; j < size; ++j){
//if(i == 0 || i == size-1
// || j == 0 || j == size-1)
// temp
//else
temp[i][j] = (i+1)*10 + j + 1;
cout << temp[i][j] << " ";
}
cout << endl;
}
cout << endl;
dct2(temp, size, size, -1);
for(int i = 0; i < size; ++i){
for(int j = 0; j < size; ++j){
cout <<setw(10)<< temp[i][j] << " ";
}
cout << endl;
}
cout << endl;
dct2(temp, size, size, 1);
for(int i = 0; i < size; ++i){
for(int j = 0; j < size; ++j){
cout <<setw(10)<< temp[i][j] << " ";
}
cout << endl;
}
cout << endl;
return 0;
}
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page