ophKernel.cu

Go to the documentation of this file.

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install, copy or use the software.
//
//
//                           License Agreement
//                For Open Source Digital Holographic Library
//
// Openholo library is free software;
// you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
//
// Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
// E-mail : contact.openholo@gmail.com
// Web : http://www.openholo.org
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//  1. Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//  2. Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the copyright holder or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
// This software contains opensource software released under GNU Generic Public License,
// NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
// Check whether software you use contains licensed software.
//
//M*/
#include "ophKernel.cuh"

__global__ void fftShift(int N, int nx, int ny, cufftDoubleComplex* input, cufftDoubleComplex* output, bool bNormalized)
{
   int tid = threadIdx.x + blockIdx.x * blockDim.x;

   double normalF = 1.0;
   if (bNormalized == true)
       normalF = nx * ny;

   int half_nx = nx >> 1;
   int half_ny = ny >> 1;
   unsigned int size = blockDim.x * gridDim.x;

   while (tid < N)
   {
       int i = tid % nx;
       int j = tid / nx;

       int ti = i - half_nx; if (ti < 0) ti += nx;
       int tj = j - half_ny; if (tj < 0) tj += ny;

       int oindex = tj * nx + ti;


       output[tid].x = input[oindex].x / normalF;
       output[tid].y = input[oindex].y / normalF;

       tid += size;
   }
}

__global__ void fftShiftf(int N, int nx, int ny, cuFloatComplex* input, cuFloatComplex* output, bool bNormalized)
{
   int tid = threadIdx.x + blockIdx.x * blockDim.x;

   float normalF = 1.0;
   if (bNormalized == true)
       normalF = nx * ny;

   int half_nx = nx >> 1;
   int half_ny = ny >> 1;

   unsigned int size = blockDim.x * gridDim.x;

   while (tid < N)
   {
       int i = tid % nx;
       int j = tid / nx;

       int ti = i - half_nx; if (ti < 0) ti += nx;
       int tj = j - half_ny; if (tj < 0) tj += ny;

       int oindex = tj * nx + ti;


       output[tid].x = input[oindex].x / normalF;
       output[tid].y = input[oindex].y / normalF;

       tid += size;
   }
}

__device__ void exponent_complex(cuDoubleComplex* val)
{
   double exp_val = exp(val->x);
   //double re = val->x;
   double im = val->y;

   val->x = exp_val * cos(im);
   val->y = exp_val * sin(im);
}

extern "C"
{
   void cudaFFT(CUstream_st* stream, int nx, int ny, cufftDoubleComplex* in_field, cufftDoubleComplex* output_field, int direction, bool bNormalized)
   {
       unsigned int nblocks = (nx * ny + kBlockThreads - 1) / kBlockThreads;
       int N = nx * ny;
       fftShift << <nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, false);

       cufftHandle plan;
       cufftResult result;
       // fft
       result = cufftPlan2d(&plan, ny, nx, CUFFT_Z2Z);
       if (result != CUFFT_SUCCESS)
       {
           LOG("<FAILED> cufftPlan2d (%d)\n", result);
           return;
       };

       if (direction == -1)
           result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_FORWARD);
       else
           result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_INVERSE);

       if (result != CUFFT_SUCCESS)
       {
           LOG("<FAILED> cufftExecZ2Z (%d)\n", result);
           return;
       }

       if (cudaDeviceSynchronize() != cudaSuccess) {
           LOG("<FAILED> cudaDeviceSynchronize\n");
           return;
       }

       fftShift << < nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, bNormalized);

       cufftDestroy(plan);
   }
}