Получение фазового изображения из CUDA FFT - PullRequest
Новая
       84

Получение фазового изображения из CUDA FFT

0 голосов
/ 10 октября 2018

Я пытаюсь применить cuFFT, вперед, а затем обратно, к 2D-изображению.Мне нужны реальные и сложные части в качестве отдельных выходов, чтобы я мог рассчитать изображение фазы и величины.Я не смог воссоздать входное изображение, а также возвращается ненулевая фаза.В частности, я не уверен, правильно ли я создаю полноразмерное изображение из комплексного вывода cuFFT уменьшенного размера, которое, по-видимому, хранит только левую сторону спектра.Вот мой текущий код: 

// Load image
cv::Mat_<float> img;
img = cv::imread(path,0);
if(!img.isContinuous()){
    std::cout<<"Input cv::Mat is not continuous!"<<std::endl;
    return -1;
}

float *h_Data, *d_Data;
h_Data = img.ptr<float>(0);

// Complex device pointers
cufftComplex
*d_DataSpectrum,
*d_Result,
*h_Result;

// Plans for cuFFT execution
cufftHandle
fftPlanFwd,
fftPlanInv;

// Image dimensions
const int dataH = img.rows;
const int dataW = img.cols;
const int complexW = dataW/2+1;

// Allocate memory
h_Result = (cufftComplex *)malloc(dataH    * complexW * sizeof(cufftComplex));
checkCudaErrors(cudaMalloc((void **)&d_DataSpectrum,   dataH * complexW * sizeof(cufftComplex)));
checkCudaErrors(cudaMalloc((void **)&d_Data,   dataH   * dataW   * sizeof(float)));
checkCudaErrors(cudaMalloc((void **)&d_Result,   dataH * complexW * sizeof(cufftComplex)));

// Copy image to GPU
checkCudaErrors(cudaMemcpy(d_Data,   h_Data,   dataH   * dataW *   sizeof(float), cudaMemcpyHostToDevice));

// Forward FFT
checkCudaErrors(cufftPlan2d(&fftPlanFwd, dataH, dataW, CUFFT_R2C));
checkCudaErrors(cufftExecR2C(fftPlanFwd, (cufftReal *)d_Data, (cufftComplex *)d_DataSpectrum));

// Inverse FFT
checkCudaErrors(cufftPlan2d(&fftPlanInv, dataH, dataW, CUFFT_C2C));
checkCudaErrors(cufftExecC2C(fftPlanInv, (cufftComplex *)d_DataSpectrum, (cufftComplex *)d_Result, CUFFT_INVERSE));

// Copy result to host memory
checkCudaErrors(cudaMemcpy(h_Result, d_Result, dataH * complexW * sizeof(cufftComplex), cudaMemcpyDeviceToHost));

// Convert cufftComplex to OpenCV real and imag Mat
Mat_<float> resultReal = Mat_<float>(dataH, dataW);
Mat_<float> resultImag = Mat_<float>(dataH, dataW);
for(int i=0; i<dataH; i++){
    float* rowPtrReal = resultReal.ptr<float>(i);
    float* rowPtrImag = resultImag.ptr<float>(i);
    for(int j=0; j<dataW; j++){
        if(j<complexW){
            rowPtrReal[j] = h_Result[i*complexW+j].x/(dataH*dataW);
            rowPtrImag[j] = h_Result[i*complexW+j].y/(dataH*dataW);
        }else{
            // Right side?
            rowPtrReal[j] = h_Result[i*complexW+(dataW-j)].x/(dataH*dataW);
            rowPtrImag[j] = -h_Result[i*complexW+(dataW-j)].y/(dataH*dataW);
        }
    }
}

// Compute phase and normalize to 8 bit
Mat_<float> resultPhase;
phase(resultReal, resultImag, resultPhase);
cv::subtract(resultPhase, 2*M_PI, resultPhase, (resultPhase > M_PI));
resultPhase = ((resultPhase+M_PI)*255)/(2*M_PI);
Mat_<uchar> normalized = Mat_<uchar>(dataH, dataW);
resultPhase.convertTo(normalized, CV_8U);
// Save phase image
cv::imwrite("cuda_propagation_phase.png",normalized);

// Compute amplitude and normalize to 8 bit
Mat_<float> resultAmplitude;
magnitude(resultReal, resultImag, resultAmplitude);
Mat_<uchar> normalizedAmplitude = Mat_<uchar>(dataH, dataW);
resultAmplitude.convertTo(normalizedAmplitude, CV_8U);
// Save phase image
cv::imwrite("cuda_propagation_amplitude.png",normalizedAmplitude);

Я не уверен, где моя ошибка.Это правильный способ вернуть все изображение из уменьшенной версии (цикл for)?

Ответы [ 2 ]

0 голосов
/ 25 января 2019

Это старый вопрос, но я хотел бы предоставить дополнительную информацию: R2C сохраняет тот же объем информации, что и преобразование C2C, он просто делает это примерно с половиной элементов.Преобразования R2C (и C2R) используют эрмитову симметрию, чтобы уменьшить количество элементов, которые вычисляются и хранятся в памяти (например, БПФ является симметричным, поэтому вам фактически не нужна половина терминов, хранящихся вПреобразование C2C).

Чтобы создать двумерное изображение действительных и мнимых компонентов, можно использовать преобразование R2C, а затем написать ядро, которое переводит выходной массив (Nx / 2 + 1) Ny впара массивов размера (Nx Ny), используя симметрию самостоятельно, чтобы записать термины в правильные позиции.Но использование преобразования C2C немного меньше кода и более надежно.

0 голосов
/ 11 октября 2018

Я думаю, что получил это сейчас.«Хитрость» - начать со сложной матрицы.Начиная с реального, вам нужно применить преобразование R2C - которое использует уменьшенный размер из-за симметрии спектра - и затем преобразование C2C, которое сохраняет этот уменьшенный размер.Решение состояло в том, чтобы создать сложный входной сигнал из реального, вставив нули как сложную часть, а затем применив два преобразования C2C в ряд, которые одновременно сохраняют все изображение и облегчают последующее получение полноразмерных действительных и мнимых матриц: 

// Load image
cv::Mat_<float> img;
img = cv::imread(path,0);
if(!img.isContinuous()){
    std::cout<<"Input cv::Mat is not continuous!"<<std::endl;
    return -1;
}

float *h_DataReal = img.ptr<float>(0);
cufftComplex *h_DataComplex;

// Image dimensions
const int dataH = img.rows;
const int dataW = img.cols;

// Convert real input to complex
h_DataComplex = (cufftComplex *)malloc(dataH    * dataW * sizeof(cufftComplex));
for(int i=0; i<dataH*dataW; i++){
    h_DataComplex[i].x = h_DataReal[i];
    h_DataComplex[i].y = 0.0f;
}

// Complex device pointers
cufftComplex
*d_Data,
*d_DataSpectrum,
*d_Result,
*h_Result;

// Plans for cuFFT execution
cufftHandle
fftPlanFwd,
fftPlanInv;

// Allocate memory
h_Result = (cufftComplex *)malloc(dataH    * dataW * sizeof(cufftComplex));
checkCudaErrors(cudaMalloc((void **)&d_DataSpectrum,   dataH * dataW * sizeof(cufftComplex)));
checkCudaErrors(cudaMalloc((void **)&d_Data,   dataH   * dataW   * sizeof(cufftComplex)));
checkCudaErrors(cudaMalloc((void **)&d_Result,   dataH * dataW * sizeof(cufftComplex)));

// Copy image to GPU
checkCudaErrors(cudaMemcpy(d_Data,   h_DataComplex,   dataH   * dataW *   sizeof(cufftComplex), cudaMemcpyHostToDevice));

// Forward FFT
checkCudaErrors(cufftPlan2d(&fftPlanFwd, dataH, dataW, CUFFT_C2C));
checkCudaErrors(cufftExecC2C(fftPlanFwd, (cufftComplex *)d_Data, (cufftComplex *)d_DataSpectrum, CUFFT_FORWARD));

// Inverse FFT
checkCudaErrors(cufftPlan2d(&fftPlanInv, dataH, dataW, CUFFT_C2C));
checkCudaErrors(cufftExecC2C(fftPlanInv, (cufftComplex *)d_DataSpectrum, (cufftComplex *)d_Result, CUFFT_INVERSE));

// Copy result to host memory
checkCudaErrors(cudaMemcpy(h_Result, d_Result, dataH * dataW * sizeof(cufftComplex), cudaMemcpyDeviceToHost));

// Convert cufftComplex to OpenCV real and imag Mat
Mat_<float> resultReal = Mat_<float>(dataH, dataW);
Mat_<float> resultImag = Mat_<float>(dataH, dataW);
for(int i=0; i<dataH; i++){
    float* rowPtrReal = resultReal.ptr<float>(i);
    float* rowPtrImag = resultImag.ptr<float>(i);
    for(int j=0; j<dataW; j++){
            rowPtrReal[j] = h_Result[i*dataW+j].x/(dataH*dataW);
            rowPtrImag[j] = h_Result[i*dataW+j].y/(dataH*dataW);
    }
}
...