Почему эта реализация Foward DCT версии Лоффлера не работает?
Я пытаюсь реализовать версию 1D DCT Лоффлера, но безрезультатно … Я следовал цепочке операций, показанной на блок-схеме, но изображение становится белым 🙁 Что я делаю не так?
Диаграмма:
Код:
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <stdio.h>
#include <opencv/cv.hpp>
#include <pthread.h>
#include "highgui.h"#include "cv.h"
using namespace std;
using namespace cv;
#define C3 851
#define S3 569
#define C1 1004
#define S1 200
#define C6 392
#define S6 946
#define R2 181
void dct2(Mat in, double dct[8][8], int xpos, int ypos) {
int i;
double rows[8][8];
int x0, x1, x2, x3, x4, x5, x6, x7, x8;
//cout << S3 << " " << S1 << " " << S6 << endl;
/* transform rows */
for (i = 0; i < 8; i++) {
x0 = in.at<uchar>(xpos + 0, ypos + i);
x1 = in.at<uchar>(xpos + 1, ypos + i);
x2 = in.at<uchar>(xpos + 2, ypos + i);
x3 = in.at<uchar>(xpos + 3, ypos + i);
x4 = in.at<uchar>(xpos + 4, ypos + i);
x5 = in.at<uchar>(xpos + 5, ypos + i);
x6 = in.at<uchar>(xpos + 6, ypos + i);
x7 = in.at<uchar>(xpos + 7, ypos + i);
//STAGE 1
int X0 = x0;
x0 += x7;
int X1 = x1;
x1 += x6;
int X2 = x2;
x2 += x5;
int X3 = x3;
x3 += x4;
x4 = X3 - x4;
x5 = X2 - x5;
x6 = X1 - x6;
x7 = X0 - x7;
//STAGE 2
X0 = x0;
X1 = x1;
x0 += x3;
x1 += x2;
x2 = X1 - x2;
x3 = X0 - x3;
int X4 = x4;
x4 = x4 * C3 + x7 * S3;
x7 = x7 * C3 - X4 * S3;
int X5 = x5;
x5 = x5 * C1 + x6 * S1;
x6 = x6 * C1 - X5 * S1;
//STAGE 3
X0 = x0;
x0 += x1;
x1 = X0 - x1;
X2 = x2;
x2 = R2 * (x2 * C6 + x3 * S6);
x3 = R2 * (x3 * C6 - X2 * S6);
X4 = x4;
X5 = x5;
x4 += x6;
x5 = x7 - x5;
x6 = X4 - x6;
x7 += X5;
//STAGE 4
X4 = x4;
rows[i][0] = x0;
rows[i][4] = x1;
rows[i][2] = x2 >> 17;
rows[i][6] = x3 >> 17;
rows[i][7] = (x4 + x7) >> 10;
rows[i][3] = (x5 * R2) >> 17;
rows[i][5] = (x6 * R2) >> 17;
rows[i][2] = (x4 - x7) >> 10;
}
/* transform columns */
for (i = 0; i < 8; i++) {
x0 = rows[0][i];
x1 = rows[1][i];
x2 = rows[2][i];
x3 = rows[3][i];
x4 = rows[4][i];
x5 = rows[5][i];
x6 = rows[6][i];
x7 = rows[7][i];
//STAGE 1
int X0 = x0;
x0 += x7;
int X1 = x1;
x1 += x6;
int X2 = x2;
x2 += x5;
int X3 = x3;
x3 += x4;
x4 = X3 - x4;
x5 = X2 - x5;
x6 = X1 - x6;
x7 = X0 - x7;
//STAGE 2
X0 = x0;
X1 = x1;
x0 += x3;
x1 += x2;
x2 = X1 - x2;
x3 = X0 - x3;
int X4 = x4;
x4 = x4 * C3 + x7 * S3;
x7 = x7 * C3 - X4 * S3;
int X5 = x5;
x5 = x5 * C1 + x6 * S1;
x6 = x6 * C1 - X5 * S1;
//STAGE 3
X0 = x0;
x0 += x1;
x1 = X0 - x1;
X2 = x2;
x2 = R2 * (x2 * C6 + x3 * S6);
x3 = R2 * (x3 * C6 - X2 * S6);
X4 = x4;
X5 = x5;
x4 += x6;
x5 = x7 - x5;
x6 = X4 - x6;
x7 += X5;
//STAGE 4
X4 = x4;
dct[0][i] = x0;
dct[4][i] = x1;
dct[2][i] = x2 >> 17;
dct[6][i] = x3 >> 17;
dct[7][i] = (x4 + x7) >> 10;
dct[3][i] = (x5 * R2) >> 17;
dct[5][i] = (x6 * R2) >> 17;
dct[1][i] = (x4 - x7) >> 10;
}
}
#define COEFFS(Cu,Cv,u,v) { \
if (u == 0) Cu = 1.0 / sqrt(2.0); else Cu = 1.0; \
if (v == 0) Cv = 1.0 / sqrt(2.0); else Cv = 1.0; \
}
void idct2(Mat in, double data[8][8], const int xpos, const int ypos) {
int u, v, x, y;
/* iDCT */
for (y = 0; y < 8; y++)
for (x = 0; x < 8; x++) {
double z = 0.0;
for (v = 0; v < 8; v++)
for (u = 0; u < 8; u++) {
double S, q;
double Cu, Cv;
COEFFS(Cu, Cv, u, v);
S = data[v][u];
q = Cu * Cv * S
* cos(
(double) (2 * x + 1) * (double) u * M_PI
/ 16.0)
* cos(
(double) (2 * y + 1) * (double) v * M_PI
/ 16.0);
z += q;
}
z /= 4.0;
if (z > 255.0)
z = 255.0;
if (z < 0)
z = 0.0;
in.at<uchar>(x + xpos, y + ypos) = (uchar) z;
}
}
int main() {
Mat in = imread("lena.bmp", CV_LOAD_IMAGE_GRAYSCALE);
double DCT[8][8];
for (int x = 0; x < 8; ++x) {
for (int y = 0; y < 8; ++y) {
dct2(in, DCT, x * 8, y * 8);
idct2(in, DCT, x * 8, y * 8);
}
}
imshow("original", in);
waitKey(0);
return 0;
}
Результат фото Лены 64 х 64 пикселя:
Решение
Я получил! Никто не помогает мне. Просто идиот, который называет меня ленивым. Но я не! Вот код Может помочь кому-то …
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <stdio.h>
#include <opencv/cv.hpp>
#include <pthread.h>
#include "highgui.h"#include "cv.h"
using namespace std;
using namespace cv;
/*
#define C1 cos(M_PI/16)
#define C3 cos(3*M_PI/16)
#define C5 cos(5*M_PI/16)
#define C6 cos(6*M_PI/16)
#define S6 sin(6*M_PI/16)
#define C7 cos(7*M_PI/16)
#define R2 sqrt(2.0)
*/
#define C1 1004
#define C3 851
#define C5 569
#define C6 392
#define S6 946
#define C7 200
#define R2 181
void dct2(Mat in, double dct[8][8], int xpos, int ypos) {
int i;
double rows[8][8];
int x0, x1, x2, x3, x4, x5, x6, x7;
int tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
int tmp10, tmp11, tmp12, tmp13;
int z1, z2, z3, z4, z5;
/* transform rows */
for (i = 0; i < 8; i++) {
x0 = in.at<uchar>(xpos + 0, ypos + i);
x1 = in.at<uchar>(xpos + 1, ypos + i);
x2 = in.at<uchar>(xpos + 2, ypos + i);
x3 = in.at<uchar>(xpos + 3, ypos + i);
x4 = in.at<uchar>(xpos + 4, ypos + i);
x5 = in.at<uchar>(xpos + 5, ypos + i);
x6 = in.at<uchar>(xpos + 6, ypos + i);
x7 = in.at<uchar>(xpos + 7, ypos + i);
tmp0 = x0 + x7;
tmp7 = x0 - x7;
tmp1 = x1 + x6;
tmp6 = x1 - x6;
tmp2 = x2 + x5;
tmp5 = x2 - x5;
tmp3 = x3 + x4;
tmp4 = x3 - x4;
tmp10 = tmp0 + tmp3;
tmp13 = tmp0 - tmp3;
tmp11 = tmp1 + tmp2;
tmp12 = tmp1 - tmp2;
rows[i][0] = tmp10 + tmp11;
rows[i][4] = tmp10 - tmp11;
rows[i][2] = (R2 * (tmp12 * C6 + tmp13 * S6)) >> 17;
rows[i][6] = (R2 * (tmp13 * C6 - tmp12 * S6)) >> 17;
//ODD PART
z1 = tmp4 + tmp7;
z2 = tmp5 + tmp6;
z3 = tmp4 + tmp6;
z4 = tmp5 + tmp7;
z5 = (z3 + z4) * R2 * C3;
tmp4 = tmp4 * R2 * (-C1 + C3 + C5 - C7);
tmp5 = tmp5 * R2 * (C1 + C3 - C5 + C7);
tmp6 = tmp6 * R2 * (C1 + C3 + C5 - C7);
tmp7 = tmp7 * R2 * (C1 + C3 - C5 - C7);
z1 = z1 * R2 * (C7 - C3);
z2 = z2 * R2 * (-C1 - C3);
z3 = z3 * R2 * (-C3 - C5);
z4 = z4 * R2 * (C5 - C3);
z3 += z5;
z4 += z5;
rows[i][7] = (tmp4 + z1 + z3) >> 17;
rows[i][5] = (tmp5 + z2 + z4) >> 17;
rows[i][3] = (tmp6 + z2 + z3) >> 17;
rows[i][1] = (tmp7 + z1 + z4) >> 17;
//cout << trunc(rows[i][2]) << endl;
}
/* transform columns */
for (i = 0; i < 8; i++) {
x0 = rows[0][i];
x1 = rows[1][i];
x2 = rows[2][i];
x3 = rows[3][i];
x4 = rows[4][i];
x5 = rows[5][i];
x6 = rows[6][i];
x7 = rows[7][i];
tmp0 = x0 + x7;
tmp7 = x0 - x7;
tmp1 = x1 + x6;
tmp6 = x1 - x6;
tmp2 = x2 + x5;
tmp5 = x2 - x5;
tmp3 = x3 + x4;
tmp4 = x3 - x4;
tmp10 = tmp0 + tmp3;
tmp13 = tmp0 - tmp3;
tmp11 = tmp1 + tmp2;
tmp12 = tmp1 - tmp2;
dct[0][i] = (tmp10 + tmp11) >> 3;
dct[4][i] = (tmp10 - tmp11) >> 3;
dct[2][i] = (R2 * (tmp12 * C6 + tmp13 * S6)) >> 20;
dct[6][i] = (R2 * (tmp13 * C6 - tmp12 * S6)) >> 20;
//cout << dct[0][i] << endl;
//ODD PART
z1 = tmp4 + tmp7;
z2 = tmp5 + tmp6;
z3 = tmp4 + tmp6;
z4 = tmp5 + tmp7;
z5 = (z3 + z4) * R2 * C3;
tmp4 = tmp4 * R2 * (-C1 + C3 + C5 - C7);
tmp5 = tmp5 * R2 * (C1 + C3 - C5 + C7);
tmp6 = tmp6 * R2 * (C1 + C3 + C5 - C7);
tmp7 = tmp7 * R2 * (C1 + C3 - C5 - C7);
z1 = z1 * R2 * (C7 - C3);
z2 = z2 * R2 * (-C1 - C3);
z3 = z3 * R2 * (-C3 - C5);
z4 = z4 * R2 * (C5 - C3);
z3 += z5;
z4 += z5;
dct[7][i] = (tmp4 + z1 + z3) >> 20;
dct[5][i] = (tmp5 + z2 + z4) >> 20;
dct[3][i] = (tmp6 + z2 + z3) >> 20;
dct[1][i] = (tmp7 + z1 + z4) >> 20;
//cout << dct[0][i] << endl;
}
}
#define COEFFS(Cu,Cv,u,v) { \
if (u == 0) Cu = 1.0 / sqrt(2.0); else Cu = 1.0; \
if (v == 0) Cv = 1.0 / sqrt(2.0); else Cv = 1.0; \
}
void idct2(Mat in, double data[8][8], const int xpos, const int ypos) {
int u, v, x, y;
/* iDCT */
for (y = 0; y < 8; y++)
for (x = 0; x < 8; x++) {
double z = 0.0;
for (v = 0; v < 8; v++)
for (u = 0; u < 8; u++) {
double S, q;
double Cu, Cv;
COEFFS(Cu, Cv, u, v);
S = data[v][u];
q = Cu * Cv * S
* cos(
(double) (2 * x + 1) * (double) u * M_PI
/ 16.0)
* cos(
(double) (2 * y + 1) * (double) v * M_PI
/ 16.0);
z += q;
}
z /= 4.0;
if (z > 255.0)
z = 255.0;
if (z < 0)
z = 0.0;
in.at<uchar>(x + xpos, y + ypos) = (uchar) z;
}
}
int main() {
Mat in = imread("lena.bmp", CV_LOAD_IMAGE_GRAYSCALE);
double DCT[8][8];
for (int x = 0; x < in.cols/8; ++x) {
for (int y = 0; y < in.rows/8; ++y) {
dct2(in, DCT, x * 8, y * 8);
idct2(in, DCT, x * 8, y * 8);
}
}
imshow("original", in);
waitKey(0);
return 0;
}
Другие решения
Других решений пока нет …