Список использованных источников. Какие же можно сделать выводы.

Выводы

Какие же можно сделать выводы? Первое, что хочется отметить – это явное отставание с помощью фильтрации с применением преобразования Фурье.

Вопреки ожиданиям, этот метод показал себя крайне неэффективно и заочно выбыл из испытаний. Как такое могло случиться? Возможно, это связано с огромным количеством циклических операций сложения в реализации метода. Так же может повлиять постоянная адресация в различные сегменты памяти ОЗУ, где влияет уже латентность(memory latency) самого физического устройства. Возможно, из-за этой задержки сказывается такое поведение.

Так же в этой реализации присутствуют множество операций с плавающей точкой. Хотя современные процессоры научились хорошо работать в этом режиме – но в совокупности факторов это могло повлиять. Я так же проводил тесты, задавая указание компилятору, внутренне представлять числа с плавающей точкой в виде фиксированной. Но данный метод не сильно оправдал себя. По результатам исследования я получил не более 6% прироста производительности (~6400 миллисекунд против ~6000). В качестве вывода для этого метода можно предположить следующее: без соответствующей хорошей реализации самого преобразования Фурье мы будем худшую производительность за счёт огромного количества проходов в циклах, где осуществляются операции с плавающей точкой и неструктурные запросы в оперативную память.

Явный метод показал себя, как и следовало ожидать, не в первых местах. Однако в данном исследовании не на последнем месте. Здесь сказывается большое количество операций - M*N*D * D, где М – изображение в ширину, N – изображениев высоту, D – размер окна. Как видим, здесь каких либо оптимизаций сделать не удастся, разве что за счёт уменьшения размера окна. Но вопрос, зачем, если это уменьшает эффективность самого фильтра.

Намного лучше предыдущих показал себя рекурсивный фильтр Гаусса. Здесь мы видим скорость и эффективность недоступную вышеописанным претендентам. Хотя он и имеет сложность M*N*4, за счёт прохождения по строкам изображения, по столбцам и дважды транспонирования, но показал сея как вполне эффективный фильтр, который уже можно применять в прикладных программах. Что можно улучшить в данном фильтре – так это избавиться от транспонирования, что даст небольшой, но всё же прирост производительности.

И наконец, наиболее эффективный алгоритм нашего исследования – рекурсивный фильтр Гаусса с использованием параллелизации на GPU и CPU средствами OpenCL. Данный метод показал себя лучше всех даже на тестовом двуядерном CPU. В некоторых тестах его производительность была более чем в три раза выше ближайшего конкурента. А это очень много. Как можно улучшить данный метод? Со стороны алгоритма, как и в рекурсивном случае, избавиться от транспонирования, а со стороны аппаратного обеспечения –повысить количество рабочих ядер. OpenCL позволяет задействовать все имеющиеся в наличии ресурсы компьютера, как GPU, так и CPU, поэтому, как говорится, чем больше – тем лучше. И данный метод уже может применяться на больших объёмах данных с использованием мощных кластеров.

ЗАКЛЮЧЕНИЕ

1. В работе реализованы основные алгоритмы вычисления гауссова фильтра, применяемого для сглаживания растровых изображений: явный метод, метод быстрого преобразования Фурье, рекурсивный метод. Проведен сравнительный анализ этих методов.

2. Рекурсивный метод является наиболее пригодным для эффективного сглаживания растровых изображений, не смотря на то, что он вычисляет результирующее изображение с некоторой погрешностью (0.69%). Такая погрешность допустима в подавляющем большинстве практических задач современной компьютерной графики.

3. Реализована параллелизация рекурсивного алгоритма с использованием библиотеки OpenCL. В результате быстродействие этого метода повышено в 2.5-3 раза.

1. Young I., van Vliet L. Recursive implementation of the Gaussian filter. Signal Processing, vol. 44. ­– Elsevier, 1995. – C.139-151.

2. Witkin A. Scale-space filtering. – Karlsruhe, Germany, 1983. – C. 1019-1021.

3. Shapiro L., Stockman G. Computer Vision. — Prentence Hall, 2001. — C. 137-150.

4. Nixon M., Aguado A. Feature Extraction and Image Processing. — Academic Press, 2008. — C. 88.

5. Сергиенко А. Б. Цифровая обработка сигналов. — Спб: Питер, 2006. — С. 751.

6. Papoulis A. Signal Analysis. – McGraw-Hill, 1977. – C. 23.

7. Рудин У. Основы математического анализа. – М., 1976.

ПРИЛОЖЕНИЕ А. Руководство пользователя

Для корректной работы приложения необходимо установить OpenCL на компьютер. После запуска программы необходимо открыть файл необходимо открыть какой-либо файл с высотой и шириной степени двойки:

После в меню выбрать какой либо фильтр и нажать применить:

Так же можно посмотреть один из этапов преобразования Фурье выбрав в меню соответствующий пункт:

ПРИЛОЖЕНИЕ Б. Листинг программы

#include <complex>

#include <algorithm>

#include <qmath.h>

#include <QtDebug>

#include <QTime>

#include <QImageReader>

#include "gauss.h"

#include "fourier.h"

#include "..\..\qtopencl\qclcontext.h"

#define idxc(x,y, stride) ((y)*(stride)*3+(x)*3 + c)

#define idx(x,y, stride) ((y)*(stride)+(x))

Gauss::Gauss():width(-1),height(-1),bpp(-1) {

static double predefinedMatrix[] = {

0.00000067,0.00002292,0.00019117,0.00038771,0.00019117,0.00002292,0.00000067,

0.00002292,0.00078633,0.00655965,0.01330373,0.00655965,0.00078633,0.00002292,

0.00019117,0.00655965,0.05472157,0.11098164,0.05472157,0.00655965,0.00019117,

0.00038771,0.01330373,0.11098164,0.22508352,0.11098164,0.01330373,0.00038771,

0.00019117,0.00655965,0.05472157,0.11098164,0.05472157,0.00655965,0.00019117,

0.00002292,0.00078633,0.00655965,0.01330373,0.00655965,0.00078633,0.00002292,

0.00000067,0.00002292,0.00019117,0.00038771,0.00019117,0.00002292,0.00000067

};

//we regenerate it by ourselves

genereateNativeGaussian(7, 0.84089642);

//generateFourierMatrix(5);

qDebug()<<"------------------------";

}

QVector<double> Gauss::genFou(int dim){

QVector<double> matrix(dim*dim);

double sigma = native.sigma;

for (int i = 0; i < dim; i++){

for (int j = 0; j < dim; j++){

double p = (-i * i - j * j)/(2.0 * sigma * sigma);

double e = exp(p);

matrix[ j*dim + i ] = e;

//qDebug()<<QString("%1 k=%2 p=%3 e=%4").arg(matrix[idx-1], 12, 'f',8).arg(k, 12, 'f',8).arg(p, 12, 'f',8).arg(e, 12, 'f',8);

}

}

return matrix;

}

QVector<double> Gauss::genereateNativeGaussian(int dimension, double sigma, bool f)

{

qDebug()<<"generating gaussian matrix";

if (dimension == -1)

dimension = native.dim;

if (sigma == -1)

sigma = native.sigma;

QVector<double> matrix = QVector<double>(dimension*dimension);

int half_dim = dimension/2, idx = 0;

double sum = 0;

int odd = dimension % 2;

for (int i = -half_dim; i <= half_dim - 1 + odd; i++){

for (int j = -half_dim; j <= half_dim - 1 + odd; j++){

double k = 1.0 / sqrt(2.0 * M_PI) / sigma;

double p = (-i * i - j * j)/(2.0 * sigma * sigma);

double e = exp(p);

sum += (matrix[ idx++ ] = (!f?k:1) * e);

}

}

//normalize

for (uint i = 0;!f && i < (uint)dimension * dimension; i++)

matrix[i] /= sum;

native.dim = dimension;

native.matrix = matrix;

native.sigma = sigma;

if (!"show matrix"){

showNative(native.matrix, native.dim);

QVector<double> im(dimension*dimension);

QVector<double> dRE(dimension*dimension);

QVector<double> dIM(dimension*dimension);

DFT2D<1>(dimension, dimension,0,native.matrix, im, dRE, dIM);

showNative(dRE, native.dim);

DFT2D<1>(dimension, dimension,1, dRE, dIM, native.matrix, im);

showNative(native.matrix, native.dim);

showNative(im, native.dim);

}

return native.matrix;

}

bool Gauss::openImage(QString filePath){

loaded = false;

QImageReader reader(filePath);

orig.image = reader.read();

if (!orig.image.isNull()){

loaded = true;

orig.image = orig.image.convertToFormat(QImage::Format_RGB32);

width = orig.image.width();

height = orig.image.height();

bpp = orig.image.depth();

initFourier();

} else

qDebug()<<reader.errorString();

return loaded;

}

QPixmap Gauss::originalPixmap(){

return QPixmap::fromImage(orig.image);

}

QPixmap Gauss::nativeFilter(){

return QPixmap::fromImage(native.image);

}

int Gauss::applyNativeFilter(){

QRgb *bits = (QRgb*)orig.image.bits();

int qrgbCount = orig.image.byteCount() / sizeof(QRgb);

//TODO free

QRgb *gauss = new QRgb[qrgbCount];

QTime t;

t.start();

for (int x = 0; x < width; x++){

for (int y = 0; y < height; y++){

int pixelIndex = y * width + x;

double resultRed=0;

double resultGreen=0;

double resultBlue=0;

for (int x0 = x - native.dim / 2, xg=0; x0 < x + native.dim / 2; x0++, xg++){

if (x0 < 0 || x0 >= width)

continue;

for (int y0 = y - native.dim / 2, yg=0; y0 < y + native.dim / 2; y0++, yg++){

if (y0 < 0 || y0 >= height)

continue;

double gaussian = native.matrix[ yg * native.dim + xg ];

QRgb in = bits[y0 * width + x0];

resultRed += (double)qRed(in) * gaussian;

resultGreen += (double)qGreen(in) * gaussian;

resultBlue += (double)qBlue(in) * gaussian;

}

}

gauss[pixelIndex] = qRgb(resultRed, resultGreen, resultBlue);

}

}

int elapsed = t.elapsed();

native.image = QImage((uchar*)gauss, width, height, QImage::Format_RGB32);

return elapsed;

}

int Gauss::applyFourier(int type, bool inv, bool shift)

{

QTime t;

t.start();

if (type==0)

DFT2D<3>(fourier.width2, fourier.height2, inv, fourier.fRe, fourier.fIm, fourier.FRe, fourier.FIm);

else

FFT2D<3>(fourier.width2, fourier.height2, inv, fourier.fRe, fourier.fIm, fourier.FRe, fourier.FIm);

int elapsed = t.elapsed();

fix(fourier.width2, fourier.height2, fourier.FRe, fourier.gauss, shift, false);

std::swap(fourier.fRe,fourier.FRe);

std::swap(fourier.fIm,fourier.FIm);

fourier.image = QImage((uchar*)fourier.gauss, fourier.width2, fourier.height2, QImage::Format_RGB32);

return elapsed;

}

int Gauss::applyFourierFilter()

{

QTime t;

t.start();

QVector<double> mIm;

QVector<double> mRE;

QVector<double> mIM;

if ("fft gaussian kernel"){

genereateNativeGaussian(fourier.width2, native.sigma, true);

int d = native.dim * native.dim;

mIm = QVector<double>(d);

mRE = QVector<double>(d);

mIM = QVector<double>(d);

DFT2D<1>(native.dim, native.dim, 0, native.matrix, mIm, mRE, mIM);

}

FFT2D<3>(fourier.width2, fourier.height2, 0, fourier.fRe, fourier.fIm, fourier.FRe, fourier.FIm);

std::swap(fourier.fRe,fourier.FRe);

std::swap(fourier.fIm,fourier.FIm);

if (!"do small kernel"){

QVector<double> matrix = native.matrix;

int dim = native.dim;

int cent = dim / 2;

for (int i = 0; i < dim/2; i++){

for (int j = 0; j < dim/2; j++){

for(int c = 0; c < 3; c++){

if ("fft g kernel"){

if ("do real"){

fourier.fRe[j*dim*3 + i*3 + c] *= mRE[(j)*dim + i];

fourier.fRe[j*dim*3 + (dim - i - 1)*3 + c] *= mRE[(j)*dim + (dim - i - 1)];

fourier.fRe[(dim - j - 1)*dim*3 + i*3 + c] *= mRE[(dim - j - 1)*dim + i];

fourier.fRe[(dim - j - 1)*dim*3 + (dim - i - 1)*3 + c] *= mRE[(dim - j -1)*dim + (dim - i - 1)];

}

if ("do imag"){

fourier.fIm[j*dim*3 + i*3 + c] *= mIM[(j)*dim + i ];

fourier.fIm[j*dim*3 + (dim - i - 1)*3 + c] *= mIM[(j)*dim + (dim - i - 1)];

fourier.fIm[(dim - j - 1)*dim*3 + i*3 + c] *= mIM[(dim - j - 1)*dim + i];

fourier.fIm[(dim - j - 1)*dim*3 + (dim - i - 1)*3 + c] *= mIM[(dim - j - 1)*dim + (dim - i - 1)];

}

} else {

if (!"do real"){

fourier.fRe[j*dim*3 + i*3 + c] *= matrix[(j + cent)*dim + i + cent];

fourier.fRe[j*dim*3 + (dim - i)*3 + c] *= matrix[(j + cent)*dim + (dim - i - cent)];

fourier.fRe[(dim - j - 1)*dim*3 + i*3 + c] *= matrix[(dim - j - cent)*dim + i + cent];

fourier.fRe[(dim - j - 1)*dim*3 + (dim - i)*3 + c] *= matrix[(dim - j - cent)*dim + (dim - i - cent)];

}

if (!"do imag"){

fourier.fIm[j*dim*3 + i*3 + c] *= matrix[(j + cent)*dim + i + cent];

fourier.fIm[j*dim*3 + (dim - i)*3 + c] *= matrix[(j + cent)*dim + (dim - i - cent)];

fourier.fIm[(dim - j - 1)*dim*3 + i*3 + c] *= matrix[(dim - j - cent)*dim + i + cent];

fourier.fIm[(dim - j - 1)*dim*3 + (dim - i)*3 + c] *= matrix[(dim - j - cent)*dim + (dim - i - cent)];

}

}

}

}

}

} else{

//genereateNativeGaussian(fourier.width2, native.sigma);

int half = fourier.width2/2;

QVector<double> matrix = genFou(fourier.width2);

for (int i = 0; i < fourier.height2; i++){

for (int j = 0; j < fourier.width2; j++){

for(int c = 0; c <3; c++){

fourier.fRe[ idxc(j,i,fourier.width2)] *= matrix[idx(j,i,fourier.width2)];

//fourier.fIm[ idxc(j,i,fourier.width2)] *= matrix[idx(j,i,fourier.width2)];

}

}

}

}

FFT2D<3>(fourier.width2, fourier.height2, 1, fourier.fRe, fourier.fIm, fourier.FRe, fourier.FIm);

int elapsed = t.elapsed();

fix(fourier.width2, fourier.height2, fourier.FRe,fourier.gauss, false, false);

std::swap(fourier.fRe,fourier.FRe);

std::swap(fourier.fIm,fourier.FIm);

fourier.image = QImage((uchar*)fourier.gauss,fourier.width2,fourier.height2,QImage::Format_RGB32);

return elapsed;

}

complex_t* Gauss::generateFourierMatrix(int n)

{

Q_UNUSED(n);

return 0;

}

QPixmap Gauss::fourierFilter()

{

return QPixmap::fromImage(fourier.image);

}

void Gauss::initFourier()

{

int width2 = fourier.width2 = (int)pow(2,ceil(log(double(width))/log(2.0)));

int height2 = fourier.height2 = (int)pow(2,ceil(log(double(height))/log(2.0)));

const QRgb* bits = (const QRgb*)orig.image.bits();

QRgb* gauss = fourier.gauss = new QRgb[width2*height2 ];

size_t length = width2*height2*3;

QVector<double> &fRe = fourier.fRe = QVector<double>(length);

if (width2<=2500){

fourier.fIm = QVector<double>(length);

fourier.FRe = QVector<double>(length);

fourier.FIm = QVector<double>(length);

}

size_t tt = width2*height2;

std::fill(gauss, gauss + tt, 0);

//set signal to the image

for(int x = 0; x < width2; x++)

for(int y = 0; y < height2; y++)

{

if (x>=width)

continue;

if (y>=height)

continue;

QRgb in = bits[y*width + x];

fRe[3*y*width2 + 3*x + 0] = qRed(in);

fRe[3*y*width2 + 3*x + 1] = qGreen(in);

fRe[3*y*width2 + 3*x + 2] = qBlue(in);

}

}

int Gauss::applyRecursive()

{

Q_ASSERT(native.sigma >= 0.5);

if (0 <= native.sigma && native.sigma <= 2.5)

recursive.q = 3.97156 - 4.14554 * sqrt(1.0 - 0.26891*native.sigma);

else if (native.sigma > 2.5)

recursive.q = -0.9633 + 0.98711 * native.sigma;

double q = recursive.q;

double b0 = recursive.b0 = 1.57825 + (2.44413 * q) + (1.4281 * q * q);

double b1 = recursive.b1 = (2.44413 * q) + (2.85619 * q * q) + (1.26661 * q * q * q);

double b2 = recursive.b2 = -((1.4281 * q * q) + (1.26661 * q * q * q));

double b3 = recursive.b3 = 0.422205 * q * q * q;

double B = recursive.B = 1.0 - ((b1 + b2 + b3)/b0);

QVector<double> in = fourier.fRe;

QTime t;

t.start();

recursiveStep(in.data(), B, b1, b2, b3, b0);//fixme

transpose(in.data());

recursiveStep(in.data(), B, b1, b2, b3, b0);

transpose(in.data());

int elapsed = t.elapsed();

fix(fourier.width2, fourier.height2, in,fourier.gauss, false, false);

recursive.image = QImage((uchar*)fourier.gauss,width,height,QImage::Format_RGB32);

return elapsed;

}

QPixmap Gauss::recursiveFilter()

{

return QPixmap::fromImage(recursive.image);

}

void Gauss::transpose(double * in)

{

Q_ASSERT_X(width == height,"Recursive filter","Image is not square!");

for (int y = 0; y < height; y++){

for (int x = y; x < width; x++){

std::swap(in[x*width*3 + y*3 + 0], in[y*width*3 + x*3 + 0]);

std::swap(in[x*width*3 + y*3 + 1], in[y*width*3 + x*3 + 1]);

std::swap(in[x*width*3 + y*3 + 2], in[y*width*3 + x*3 + 2]);

}

}

std::swap(width,height);

}

double * Gauss::recursiveStep(double * ind, double B, double b1, double b2, double b3, double b0)

{

for (int c = 0; c < 3; c++){

for (int y = 0; y < height; y++){

double *in = ind + y * width * 3;

double V = in[0 + c];

in[0*3 + c] = B * V + (b1 * V + b2 * V + b3*V) / b0;

in[1*3 + c] = B * in[1*3 + c] + (b1 * in[0*3 + c] + b2 * V + b3 * V) / b0;

in[2*3 + c] = B * in[2*3 + c] + (b1 * in[1*3 + c] + b2 * in[0*3 + c] + b3 * V) / b0;

for (int i = 3; i < width; ++i)

in [i*3 + c] = B * in[i*3 + c] + (b1 * in[i*3 - 1*3 + c] + b2 * in[i*3 - 2*3 + c] + b3 * in[i*3 - 3*3 + c]) / b0;

int r = width - 1;

V = in[r*3 + c];

in [r*3 - 0*3 + c] = B * V + (b1 * V + b2 * V + b3 * V) / b0;

in [r*3 - 1*3 + c] = B * in[r*3 - 1*3 + c] + (b1 * in[r*3 + c] + b2 * V + b3 * V) / b0;

in [r*3 - 2*3 + c] = B * in[r*3 - 2*3 + c] + (b1 * in[r*3 - 1*3 + c] + b2 * in[r*3 + c] + b3 * V) / b0;

for (int i = r - 3; i >= 0; --i)

in[i*3 + c] = B * in[i*3 + c] + (b1 * in[i*3 + 1*3 + c] + b2 * in[i*3 + 2*3 + c] + b3 * in[i*3 + 3*3 + c]) / b0;

}

} return 0;

}

int Gauss::applyRecursiveCL()

{

int width2 = fourier.width2, height2 = fourier.height2;

QCLContext context;

if (!context.create()){

qDebug()<<QCLContext::errorName(context.lastError());

return 0;

}

qDebug()<< context.defaultDevice().vendor()<< context.defaultDevice().version();

QCLProgram program = context.buildProgramFromSourceFile(":/kernel");

if (program.isNull()){

qDebug()<<program.log()<<QCLContext::errorName(context.lastError());

return 0;

}

QCLKernel recursiveKernel = program.createKernel("recursive_kernel");

if (recursiveKernel.isNull()){

qDebug()<<QCLContext::errorName(context.lastError());

return 0;

}

QCLKernel transposeKernel = program.createKernel("transpose");

if (transposeKernel.isNull()){

qDebug()<<QCLContext::errorName(context.lastError());

return 0;

}

//transposeKernel.setGlobalWorkSize(width2, height2);

//transposeKernel.setLocalWorkSize(8, 8);

recursiveKernel.setGlobalWorkSize(height2, 3);

Q_ASSERT(native.sigma >= 0.5);

if (0 <= native.sigma && native.sigma <= 2.5)

recursive.q = 3.97156 - 4.14554 * sqrt(1.0 - 0.26891*native.sigma);

else if (native.sigma > 2.5)

recursive.q = -0.9633 + 0.98711 * native.sigma;

double q = recursive.q;

float b0 = recursive.b0 = 1.57825 + (2.44413 * q) + (1.4281 * q * q);

float b1 = recursive.b1 = (2.44413 * q) + (2.85619 * q * q) + (1.26661 * q * q * q);

float b2 = recursive.b2 = -((1.4281 * q * q) + (1.26661 * q * q * q));

float b3 = recursive.b3 = 0.422205 * q * q * q;

float B = recursive.B = 1.0 - ((b1 + b2 + b3)/b0);

QCLVector<float> vec = context.createVector<float>(width2*height2*3);

for (int i = 0; i < width2 * height2 * 3; i++)

vec[i] = fourier.fRe[i];

QTime t;

t.start();

vec.map();

recursiveKernel(vec, B, b1, b2, b3, b0, width2, height2);

transposeKernel(vec, width2, height2);

recursiveKernel(vec, B, b1, b2, b3, b0, width2, height2);

QCLEvent clevent = transposeKernel(vec);

clevent.waitForFinished();

vec.unmap();

int elapsed = t.elapsed();

QVector<float> in = QVector<float>(width2*height2*3);

vec.read(in.data(),width2*height2*3);

fix(fourier.width2, fourier.height2, in,fourier.gauss, false, false);

recursive.image = QImage((uchar*)fourier.gauss,width2,height2,QImage::Format_RGB32);

return elapsed;

}

void Gauss::showNative(QVector<double> mx, int dim)

{

QString s = "--------------------\n";

for (int i = 0; i < dim; i++){

for (int j = 0; j < dim; j++)

s.append(QString("%1 ").arg(mx[dim * i + j], 12, 'f',8));

s += "\n";

}

qDebug()<<s;

}

#ifndef __FFT_H__

#define __FFT_H__

#include <qmath.h>

#include <QRgb>

#ifdef DEBUG

#include <QtDebug>

#endif

#include "utils.h"

int prepare(int n, int m);

template<typename int D, typename C1, typename C2, typename C3, typename C4>

void FFT2D(int n, int m, bool inverse, C1& gRe, C2& gIm, C3& GRe, C4& GIm)

{

int l2n = 0, p = 1; //l2n will become log_2(n)

while(p < n) {p *= 2; l2n++;}

int l2m = 0; p = 1; //l2m will become log_2(m)

while(p < m) {p *= 2; l2m++;}

m= 1<<l2m; n= 1<<l2n; //Make sure m and n will be powers of 2, otherwise you'll get in an infinite loop

//Erase all history of this array

for(int x = 0; x <m; x++){ //for each column

for(int y = 0; y < m; y++){ //for each row

for(int c = 0; c < D; c++) //for each color component

{

GRe[D * m * x + D * y + c] = gRe[D * m * x + D * y + c];

GIm[D * m * x + D * y + c] = gIm[D * m * x + D * y + c];

}

}

}

//Bit reversal of each row

int j;

for(int y = 0; y < m; y++){ //for each row

for(int c = 0; c < D; c++) //for each color component

{

j = 0;

for(int i = 0; i < n - 1; i++)

{

GRe[D * m * i + D * y + c] = gRe[D * m * j + D * y + c];

GIm[D * m * i + D * y + c] = gIm[D * m * j + D * y + c];

int k = n / 2;

while (k <= j) {j -= k; k/= 2;}

j += k;

}

}

}

//Bit reversal of each column

double tx = 0, ty = 0;

for(int x = 0; x < n; x++){ //for each column

for(int c = 0; c < D; c++) //for each color component

{

j = 0;

for(int i = 0; i < m - 1; i++)

{

if(i < j)

{

tx = GRe[D * m * x + D * i + c];

ty = GIm[D * m * x + D * i + c];

GRe[D * m * x + D * i + c] = GRe[D * m * x + D * j + c];

GIm[D * m * x + D * i + c] = GIm[D * m * x + D * j + c];

GRe[D * m * x + D * j + c] = tx;

GIm[D * m * x + D * j + c] = ty;

}

int k = m / 2;

while (k <= j) {j -= k; k/= 2;}

j += k;

}

}

}

//Calculate the FFT of the columns

for(int x = 0; x < n; x++){ //for each column

for(int c = 0; c < D; c++) //for each color component

{

//This is the 1D FFT:

double ca = -1.0;

double sa = 0.0;

int l1 = 1, l2 = 1;

for(int l=0;l<l2n;l++)

{

l1 = l2;

l2 *= 2;

double u1 = 1.0;

double u2 = 0.0;

for(int j = 0; j < l1; j++)

{

for(int i = j; i < n; i += l2)

{

int i1 = i + l1;

double t1 = u1 * GRe[D * m * x + D * i1 + c] - u2 * GIm[D * m * x + D * i1 + c];

double t2 = u1 * GIm[D * m * x + D * i1 + c] + u2 * GRe[D * m * x + D * i1 + c];

GRe[D * m * x + D * i1 + c] = GRe[D * m * x + D * i + c] - t1;

GIm[D * m * x + D * i1 + c] = GIm[D * m * x + D * i + c] - t2;

GRe[D * m * x + D * i + c] += t1;

GIm[D * m * x + D * i + c] += t2;

}

double z = u1 * ca - u2 * sa;

u2 = u1 * sa + u2 * ca;

u1 = z;

}

sa = sqrt((1.0 - ca) / 2.0);

if(!inverse) sa = -sa;

ca = sqrt((1.0 + ca) / 2.0);

}

}

}

//Calculate the FFT of the rows

for(int y = 0; y < m; y++) {//for each row

for(int c = 0; c < D; c++) //for each color component

{

//This is the 1D FFT:

double ca = -1.0;

double sa = 0.0;

int l1= 1, l2 = 1;

for(int l = 0; l < l2m; l++)

{

l1 = l2;

l2 *= 2;

double u1 = 1.0;

double u2 = 0.0;

for(int j = 0; j < l1; j++)

{

for(int i = j; i < n; i += l2)

{

int i1 = i + l1;

double t1 = u1 * GRe[D * m * i1 + D * y + c] - u2 * GIm[D * m * i1 + D * y + c];

double t2 = u1 * GIm[D * m * i1 + D * y + c] + u2 * GRe[D * m * i1 + D * y + c];

GRe[D * m * i1 + D * y + c] = GRe[D * m * i + D * y + c] - t1;

GIm[D * m * i1 + D * y + c] = GIm[D * m * i + D * y + c] - t2;

GRe[D * m * i + D * y + c] += t1;

GIm[D * m * i + D * y + c] += t2;

}

double z = u1 * ca - u2 * sa;

u2 = u1 * sa + u2 * ca;

u1 = z;

}

sa = sqrt((1.0 - ca) / 2.0);

if(!inverse) sa = -sa;

ca = sqrt((1.0 + ca) / 2.0);

}

}

}

int d;

if(inverse) d = n; else d = m;

for(int x = 0; x < n; x++) for(int y = 0; y < m; y++) for(int c = 0; c < D; c++) //for every value of the buffers

{

GRe[D * m * x + D * y + c] /= d;

GIm[D * m * x + D * y + c] /= d;

}

}

/********************************************/

double get_cos(double a);

double get_sin(double a);

template<typename int D, typename C1, typename C2, typename C3, typename C4>

void DFT2D(int n, int m, bool inverse, C1& gRe, C2& gIm, C3& GRe, C4& GIm){

std::vector<double> Gr2(m * n * D);

std::vector<double> Gi2(m * n * D); //temporary buffers

//calculate the fourier transform of the columns

for(int x = 0; x < n; x++){

for(int c = 0; c < D; c++)

{

//This is the 1D DFT:

for(int w = 0; w < m; w++)

{

Gr2[m * D * x + D * w + c] =Gi2[m * D * x + D * w + c] = 0;

for(int y = 0; y < m; y++)

{

double a= 2 * M_PI * w * y / double(m);

if(!inverse)a = -a;

double ca = get_cos(a);

double sa = get_sin(a);

Gr2[m * D * x + D * w + c] += gRe[m * D * x + D * y + c] * ca - gIm[m * D * x + D * y + c] * sa;

Gi2[m * D * x + D * w + c] += gRe[m * D * x + D * y + c] * sa + gIm[m * D * x + D * y + c] * ca;

}

}

}

}

//calculate the fourier transform of the rows

for(int y = 0; y < m; y++){

for(int c = 0; c < D; c++)

{

//This is the 1D DFT:

for(int w = 0; w < n; w++)

{

GRe[m * D * w + D * y + c] = GIm[m * D * w + D * y + c] = 0;

for(int x = 0; x < n; x++)

{

double a = 2 * M_PI * w * x / double(n);

if(!inverse)a = -a;

double ca = get_cos(a);

double sa = get_sin(a);

GRe[m * D * w + D * y + c] += Gr2[m * D * x + D * y + c] * ca - Gi2[m * D * x + D * y + c] * sa;

GIm[m * D * w + D * y + c] += Gr2[m * D * x + D * y + c] * sa + Gi2[m * D * x + D * y + c] * ca;

}

if(inverse)

{

GRe[m * D * w + D * y + c] /= n;

GIm[m * D * w + D * y + c] /= n;

}

else{

GRe[m * D * w + D * y + c] /= m;

GIm[m * D * w + D * y + c] /= m;

}

}

}

}

}

struct ColorRGB{

ColorRGB():r(0),g(0),b(0){};

int r,g,b;

};

template<typename T>

void fix(int n, int m, T g, QRgb * G, bool shift, bool neg128){

Q_UNUSED(neg128);

for(int x = 0; x < n; x++)

for(int y = 0; y < m; y++)

{

int x2 = x, y2 = y;

if(shift) {x2 = (x + n / 2) % n; y2 = (y + m / 2) % m;} //Shift corners to center

ColorRGB c;

//calculate c values out of the floating point buffer

c.r = int(g[3 * m * x2 + 3 * y2 + 0]);

c.g = int(g[3 * m * x2 + 3 * y2 + 1]);

c.b = int(g[3 * m * x2 + 3 * y2 + 2]);

//negative colors give confusing effects so set them to 0

if(c.r < 0) c.r = 0;

if(c.g < 0) c.g = 0;

if(c.b < 0) c.b = 0;

//set c components higher than 255 to 255

if(c.r > 255) c.r = 255;

if(c.g > 255) c.g = 255;

if(c.b > 255) c.b = 255;

//plot the pixel

G[(x)*m+y] = qRgb(c.r,c.g,c.b);

}

}

#endif

Дата добавления: 2015-07-02