Свертка изображения (взаимная корреляция) с использованием 2D cuFFT - результат сдвинут? - PullRequest
Новая
       48

Свертка изображения (взаимная корреляция) с использованием 2D cuFFT - результат сдвинут?

0 голосов
/ 23 апреля 2020

Я пытаюсь сделать 2D-FFT для взаимной корреляции между двумя изображениями: keypoint_d размером 5x5 и image_d размером 9x9. Один из способов сделать это - использовать библиотеку cuFFT.

Итак, вот шаги, которые я использовал для преобразования IN-PLACE C2 C::

  1. Добавьте 0 отступов к Pattern_img, чтобы иметь одинаковый размер по отношению к image_d: (9x9).
  2. Я создал свой 2D C2 C план.
  3. Выполнено прямое 2D-преобразование для каждого изображения.
  4. Умножение обоих изображений в частотной области с использованием соответствующего комплексного умножения.
  5. Выполнение IFFT результата умножения.
  6. Копирование на хост и отображение результата IFFT

Вот результат моих реальных значений:

image_d: (ключевая точка существует в том же месте) 

 0  0  1  0  0  0  0  0  0
 0  1  1  1  0  0  0  0  0
 1  1  1  1  1  0  0  0  0
 0  1  1  1  0  0  0  0  0
 0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0

keypoint_d: 

 0  0  1  0  0  
 0  1  1  1  0  
 1  1  1  1  1  
 0  1  1  1  0  
 0  0  1  0  0

keypoint_padd_d: 

 0  0  1  0  0  0  0  0  0
 0  1  1  1  0  0  0  0  0
 1  1  1  1  1  0  0  0  0
 0  1  1  1  0  0  0  0  0
 0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0

IFFT-результат комплексного умножения (нормализованный) следующий: 

 0  0  0  0  0  0  0  0  0
 0  0  0  1  1  1  0  0  0
 0  0  2  3  4  3  2  0  0
 0  1  3  6  6  6  3  1  0
 0  1  4  6  10 6  4  1  0
 0  1  3  6  6  6  3  1  0
 0  0  2  3  4  3  2  0  0
 0  0  0  1  1  1  0  0  0
 0  0  0  0  0  0  0  0  0

На данный момент я хочу понять поведение IFFT. Почему наибольшее значение моего IFFT смещено (на 2) относительно центра моего шаблона? Я ожидал, что у меня будут те же координаты, что и у центра шаблона в image_d [2][2]. Другими словами, я не понимаю, почему самое высокое значение находится в [4][4]?

Я использовал код здесь с дополнительной модификацией. Вот полный код: 

int N =  5;
int SIZE = N*N;
int N2 =  9;
int sizeFig = 9 * 9;
int Gmotifx[] =
{ 0, 0, 1, 0, 0,
  0, 1, 1, 1, 0,
  1, 1, 1, 1, 1,
  0, 1, 1, 1, 0,
  0, 0, 1, 0, 0, };
Complex *fg = new Complex[SIZE];
for (int i = 0; i < SIZE; i++) {
    if (Gmotifx[i]!=0)
    {
        fg[i].x = Gmotifx[i];
        fg[i].y = 0;
    }
    else
    {
        fg[i].x = 0;
        fg[i].y = 0;
    }
}
Complex *fig = new Complex[sizeFig];
for (int i = 0; i < sizeFig; i++) {
    fig[i].x = 0; // 
    fig[i].y = 0;
}
int x_offset = 0;//2; // For shift
int y_offset = 0;// 2;
for (int i = 0; i < SIZE; i++)
{
    if (Gmotifx[i] != 0)
    {
        int x = i % N; 
        int y = i / N; 
        fig[x + x_offset + (y + y_offset) * N2].x =  Gmotifx[i]; 
        fig[x + x_offset + (y + y_offset) * N2].y = 0; 
    }
}
// ... Display the arrays
cout << "......................BEFORE PADDING " << endl;
cout << "keypoint_d::" << endl;

for (int i = 0; i < N*N - 1; i = i + N)
{
    for (int j = 0; j < N; j++)
        cout << (int)fg[i + j].x << "  ";

    cout << endl;
}

cout << "----------------" << endl;
cout << "image_d ::" << endl;

for (int i = 0; i < N2*N2 - 1; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)fig[i + j].x << "  ";

    cout << endl;
}

int mem_size = sizeof(Complex)* SIZE;
int mem_size_fig = sizeof(Complex)* sizeFig;

// .. Copy from HOST to DEVICE
cufftComplex *d_signal;
gpuErrchk(cudaMalloc((void **)&d_signal, mem_size));
gpuErrchk(cudaMemcpy(d_signal, fg, mem_size, cudaMemcpyHostToDevice));

// .. Padded signal
cufftComplex *d_signal_pad;
gpuErrchk(cudaMalloc((void **)&d_signal_pad, mem_size_fig));
gpuErrchk(cudaMemset(d_signal_pad, 0, mem_size_fig));
Complex* test2 = new Complex[sizeFig];
gpuErrchk(cudaMemcpy(test2, d_signal_pad, mem_size_fig, cudaMemcpyDeviceToHost));

cout << "-----------------create & copy into a new buffer with equal size------------------------" << endl;

for (int i = 0; i < N2*N2 ; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)test2[i + j].x << "  ";

    cout << endl;
}
cout << "-------------------------------" << endl;

for (int i = 0; i < N; i++)
{
    gpuErrchk(cudaMemcpy(d_signal_pad + N2 * i  , fg + N * i, N * sizeof(Complex), cudaMemcpyHostToDevice));

}

Complex* test = new Complex[sizeFig];
gpuErrchk(cudaMemcpy(test, d_signal_pad, mem_size_fig, cudaMemcpyDeviceToHost));


for (int i = 0; i < N2*N2; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)test[i + j].x << "  ";

    cout << endl;
}

cout << "---------------------------------------------------------------" << endl;

Complex *result = new Complex[sizeFig];
cudaMemcpy(result, d_signal_pad, sizeof(Complex)*sizeFig, cudaMemcpyDeviceToHost);

cout << "......................AFTER PADDING " << endl;
cout << "keypoint_pad_d ::" << endl;

//for (int i = 0; i < sizeFig; i = i + 10)
for (int i = 0; i < N2*N2; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)result[i + j].x << "  ";
    cout << endl;
}
cout << "-------------------------------" << endl;

cufftComplex *d_filter_kernel;
gpuErrchk(cudaMalloc((void **)&d_filter_kernel, mem_size_fig));
gpuErrchk(cudaMemcpy(d_filter_kernel, fig, mem_size_fig, cudaMemcpyHostToDevice));

Complex *result2 = new Complex[sizeFig];
cudaMemcpy(result2, d_filter_kernel, mem_size_fig, cudaMemcpyDeviceToHost);

cout << "image_d ::" << endl;

for (int i = 0; i < N2*N2; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)result2[i + j].x << "  ";
    cout << endl;
}
cout << "---------------------------------------------------------------" << endl;

// .. Start FFT :: in-place
cufftHandle plan2;
cufftPlan2d(&plan2, N2, N2, CUFFT_C2C);


printf("Transforming signal IN-PLACE cufftExecC2C\n");
cufftExecC2C(plan2, (cufftComplex *)d_signal_pad, (cufftComplex *)d_signal_pad, CUFFT_FORWARD);
cufftExecC2C(plan2, (cufftComplex *)d_filter_kernel, (cufftComplex *)d_filter_kernel, CUFFT_FORWARD);

printf("Launching Complex multiplication<<< >>>\n");
ComplexMulAndScale << < 100, 1 >> >(d_signal_pad, d_filter_kernel, 100, 1.0f / 100);

printf("Transforming signal IN-PLACE back cufftExecC2C\n");
cufftExecC2C(plan2, (cufftComplex *)d_signal_pad, (cufftComplex *)d_signal_pad, CUFFT_INVERSE);

cudaMemcpy(result2, d_signal_pad, sizeof(Complex)*sizeFig, cudaMemcpyDeviceToHost);


cout << "RESULT FFT INVERSE" << endl;
for (int i = 0; i < sizeFig; i = i + N2)
{
    for (int j = 0; j < N2; j++)
        cout << (int)result2[i+j].x << " \t ";

    cout << endl;
}

// Find the coordinates of the max value
// FindMax(result2, sizeFig, N2, N2) ; 

delete result, result2, test, test2, fg, fig;
cufftDestroy(plan2);
cudaFree(d_signal);
cudaFree(d_signal_pad);
cudaFree(d_filter_kernel);
cout << "---------------------------------------------------------------" << endl;

Ядро для сложного мульт: 

__global__ void ComplexMulAndScale (cufftComplex *a, cufftComplex *b, int size, float scale)
{
    const int tId = blockIdx.x * blockDim.x + threadIdx.x;
    if(tId < size)
    {
        Complex c;
        c.x = (a[tId].x * b[tId].x - a[tId].y * b[tId].y)*scale;
        c.y = (a[tId].x * b[tId].y + a[tId].y * b[tId].x)*scale;
        a[tId] = c;
    }
}
...