Блог

ОТКАЗ ОТ ОТВЕТСТВЕННОСТИ: Я знаю, что могу сглаживать/размывать изображения с помощью OpenCV, но мое задание более скрыто, поэтому я могу понять процесс. Я новичок в обработке изображений и DFT в целом, так что это оказывается непростой задачей.

Описание:

Я выполняю некоторую обработку изображений с использованием библиотеки FFTW и OpenCV для считывания изображения в OpenCV::Mat объект и цветовое пространство в оттенках серого. Я преобразую данные в тип double, чтобы создать на них двойной указатель, который требуется для функции FFTW для моего FFT. Это находится в объекте Mat с именем imageDouble и указателем на него pImageDouble .

 std::string filename = "<path_to_your_image>";
Mat image = imread(filename, IMREAD_GRAYSCALE);

/***********************************************************************
*** creating Mat object of type *doube* and copying image contents to it
*** then creating a pointer to it.
***********************************************************************/
Mat imageDouble;
image.convertTo(imageDouble, CV_64FC1);
double* pImageDouble = imageDouble.ptr<double>(0);
 

У меня также есть ядро Гаусса, которое я прочитал в Mat объекте. Я выполняю круговой сдвиг, при котором центральное (и самое высокое) значение ядра Гаусса смещается в верхний левый угол (положение (0,0)) ядра, а затем я обнуляю ядро до размера исходного изображения, которое я прочитал. Результаты (тип double) хранятся в объекте Mat с именем paddedGKernel и указателем с именем pPaddedGKernel .

 Mat paddedGKernel = padGKernelWithZeros(nRows, nCols, rolledGKernel);
double* pPaddedGKernel = paddedGKernel.ptr<double>(0);
 

Я инициализирую fftw_complex объекты для вывода результатов FFT в, an fftw_plan , а затем выделяю память для объектов fftw_complex и выполняю план. Я использую fftw_plan_dft_r2c_2d() функцию FFTW для выполнения БПФ для обоих 2d-объектов (изображения и смещенного/дополненного ядра Гаусса), затем выполняю точечное умножение в пространстве Фурье, чтобы применить фильтр Гаусса к исходному изображению.

 fftw_complex* out;               //for result of FFT for original image
fftw_complex* outGaussian;       //for result of FFT for Gaussian kernel
fftw_plan p;

/*************************************************
*** Allocating memory for the fftw_complex objects
*************************************************/
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);
outGaussian = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);


/****************************************
*** FFT on imageDouble, outputting to out
****************************************/
p = fftw_plan_dft_r2c_2d(imageDouble.rows, imageDouble.cols, pImageDouble, out, FFTW_ESTIMATE);
fftw_execute(p);


/**************************************************
*** FFT on paddedGKernel, outputting to outGuassian
**************************************************/
p = fftw_plan_dft_r2c_2d(paddedGKernel.rows, paddedGKernel.cols, pPaddedGKernel, outGaussian, FFTW_ESTIMATE);
fftw_execute(p);


/****************************************************
*** Pointwise multiplication to apply Gaussian kernel
****************************************************/
for (size_t i = 0; i < nCCols * nRows; i  )
{
    out[i][0] = out[i][0] * (1 / outGaussian[i][0]);
}
 

Затем я выполняю обратный БПФ fftw_plan_dft_c2r_2d() в результате точечного умножения и выводю его в объект OpenCV::Mat imageCloneDouble и преобразую данные обратно в uchar, помещая данные OpenCV::Mat imageBack в.

 /***********************************************************************
*** Making Mat object (type *double*) to put data into and pointer to it
***********************************************************************/
Mat imageCloneDouble = Mat::zeros(nRows, nCols, CV_64FC1);
double* pImageCloneDouble = imageCloneDouble.ptr<double>(0);

/*****************************************************************************
*** Making and executing plan for inverse FFT, outputting to imageCloneDouble,
*** then normalizing since this function puts out unnormalized values.
*****************************************************************************/
fftw_plan pp = fftw_plan_dft_c2r_2d(imageCloneDouble.rows, imageCloneDouble.cols, out, pImageCloneDouble, FFTW_BACKWARD | FFTW_ESTIMATE);
fftw_execute(pp);
imageCloneDouble = imageCloneDouble / (nCols * nRows *2);

/***************************************************************************
*** Making Mat object (type *uchar*) and copying data inverse FFT data to it
***************************************************************************/
Mat imageBack;
imageCloneDouble.convertTo(imageBack, CV_8UC1);
 

Когда я покажу преобразованное изображение imageBack Я ожидаю получить оригинал с примененным ядром Гаусса, но он выглядит как оригинал с потенциально некоторой модификацией, наложенной повернутой версией исходного изображения, которая, по-видимому, имеет некоторые изменения. Я не могу понять, где и почему происходит это переворачивание/вращение и почему оно накладывается на то, что кажется исходным изображением, но я подозреваю, что это происходит, когда я нахожусь в пространстве Фурье и выполняю точечное или поэлементное умножение. Либо это, либо я опускаю процесс, который мне нужно выполнить в результате точечного умножения.

Что мне нужно сделать с этими данными в пространстве Фурье, чтобы вернуться к исходному изображению?

Это твоя проблема:

Я выполняю круговой сдвиг, при котором центральное (и самое высокое) значение ядра Гаусса смещается в верхний левый угол (положение (0,0)) ядра, а затем я обнуляю ядро до размера исходного изображения, которое я прочитал.

Вам нужно изменить эти две операции на противоположные: сначала выровнять площадку до нужного размера, затем применить круговой сдвиг. Вам действительно нужно, чтобы центр ядра находился на (0,0). Но вам нужно, чтобы ядро было непрерывным в циклически периодическом мире (т. Е. При воспроизведении изображения вы должны видеть полные ядра Гаусса). Когда вы нажимаете после кругового сдвига, вы разделяете четыре квадранта ядра в этом циклически периодическом мире.

Вторая проблема заключается в этой строке:
out[i][0] = out[i][0] * (1 / outGaussian[i][0]);

Во-первых, почему делить, а не умножать? Вы пытаетесь применить свертку, нет?

Во-вторых, вы только умножаете действительную составляющую комплексных чисел, оставляя мнимую составляющую неизменной. Вам нужно выполнить полное комплексное умножение двух комплексных чисел, чтобы получить новое комплексное число. Наведите out указатель на std::complex<double>* , а затем используйте знания компилятора о сложной арифметике, чтобы выполнить ваши требования!

 std::size_t N = nRows * nCols;
fftw_complex* out = fftw_alloc_complex(N);
fftw_complex* outGaussian = fftw_alloc_complex(N);

// ...

auto* out_p = reinterpret_cast<std::complex<double>*>(out);
auto* outGaussian_p = reinterpret_cast<std::complex<double>*>(outGaussian);
for (std::size_t i = 0; i < N; i  )
{
    out_p[i] *= outGaussian_p[i];
}
 

Обратите внимание, что std::complex<double> и fftw_complex имеют одинаковую компоновку памяти , гарантированную спецификацией C для std::complex , с намерением иметь возможность выполнять этот тип приведения. На самом деле никакой переосмысления не делается, так думает только система типов C .

Ключ к устранению этой проблемы начался с того, что Крис назвал второй проблемой. Поэтому первое, что нужно было сделать, — это реализовать код, приведенный в этой части ответа.

 auto* out_p = reinterpret_cast<std::complex<double>*>(out);
auto* outGaussian_p = reinterpret_cast<std::complex<double>*>(outGaussian);
for (std::size_t i = 0; i < N; i  )
{
    out_p[i] *= outGaussian_p[i];
}
 

Я нормализовал ядро Гаусса обычным способом (разделив на сумму всех значений). После внесения изменений, указанных выше в Cris, на изображении было применено размытие, и единственная проблема заключалась в том, что изменился диапазон значений пикселей.

Чтобы исправить это, я применил формулу для растяжения контраста в соответствии с исходным диапазоном, и это решило проблему.

 /* Finding the original image's minimum and maximum value */
double oMin, oMax;
minMaxLoc(image, amp;oMin, amp;oMax);

/* Finding min and max values for our 'imageBack' and rounding them. */
double min, max;
minMaxLoc(imageCloneDouble, amp;min, amp;max);
min = round(min);
max = round(max);

/* Code version of formula for contrast stretching */
if (oMin != min | oMax != max) {
    double numerator = oMax - oMin;
    double denominator = max - min;
    double scaledMin = -(min)*numerator;
    double scaledOMin = oMin * denominator;
    double innerTerm = scaledMin   scaledOMin;

    imageCloneDouble = ((numerator * imageCloneDouble)   innerTerm ) / denominator;
}