Openholo  v4.2
Open Source Digital Holographic Library
ophKernel.cuh
Go to the documentation of this file.
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install, copy or use the software.
7 //
8 //
9 // License Agreement
10 // For Open Source Digital Holographic Library
11 //
12 // Openholo library is free software;
13 // you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
14 //
15 // Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
16 // E-mail : contact.openholo@gmail.com
17 // Web : http://www.openholo.org
18 //
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
21 //
22 // 1. Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
24 //
25 // 2. Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other materials provided with the distribution.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the copyright holder or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 // This software contains opensource software released under GNU Generic Public License,
41 // NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
42 // Check whether software you use contains licensed software.
43 //
44 //M*/
45 #pragma once
46 #ifndef __ophKernel_cuh__
47 #define __ophKernel_cuh__
48 
49 #include <cuComplex.h>
50 #include <cufft.h>
51 #include <cuda.h>
52 #include <device_launch_parameters.h>
53 #include <cuda_runtime_api.h>
54 #include "sys.h"
55 
56 static const int kBlockThreads = 512;
57 
58 inline __global__ void fftShift(int N, int nx, int ny, cufftDoubleComplex* input, cufftDoubleComplex* output, bool bNormalized)
59 {
60  int tid = threadIdx.x + blockIdx.x * blockDim.x;
61 
62  double normalF = 1.0;
63  if (bNormalized)
64  normalF = nx * ny;
65 
66  int half_nx = nx >> 1;
67  int half_ny = ny >> 1;
68  unsigned int size = blockDim.x * gridDim.x;
69 
70  while (tid < N)
71  {
72  int i = tid % nx;
73  int j = tid / nx;
74 
75  int ti = i - half_nx; if (ti < 0) ti += nx;
76  int tj = j - half_ny; if (tj < 0) tj += ny;
77 
78  int oindex = tj * nx + ti;
79 
80 
81  output[tid].x = input[oindex].x / normalF;
82  output[tid].y = input[oindex].y / normalF;
83 
84  tid += size;
85  }
86 }
87 
88 inline __global__ void fftShiftf(int N, int nx, int ny, cuFloatComplex* input, cuFloatComplex* output, bool bNormalized)
89 {
90  int tid = threadIdx.x + blockIdx.x * blockDim.x;
91 
92  float normalF = 1.0;
93  if (bNormalized)
94  normalF = nx * ny;
95 
96  int half_nx = nx >> 1;
97  int half_ny = ny >> 1;
98 
99  unsigned int size = blockDim.x * gridDim.x;
100 
101  while (tid < N)
102  {
103  int i = tid % nx;
104  int j = tid / nx;
105 
106  int ti = i - half_nx; if (ti < 0) ti += nx;
107  int tj = j - half_ny; if (tj < 0) tj += ny;
108 
109  int oindex = tj * nx + ti;
110 
111 
112  output[tid].x = input[oindex].x / normalF;
113  output[tid].y = input[oindex].y / normalF;
114 
115  tid += size;
116  }
117 }
118 
119 inline __device__ void exponent_complex(cuDoubleComplex* val)
120 {
121  double exp_val = exp(val->x);
122  //double re = val->x;
123  double im = val->y;
124 
125  val->x = exp_val * cos(im);
126  val->y = exp_val * sin(im);
127 }
128 
129 extern "C"
130 {
131  inline void cudaFFT(CUstream_st* stream, int nx, int ny, cufftDoubleComplex* in_field, cufftDoubleComplex* output_field, int direction, bool bNormalized)
132  {
133  unsigned int nblocks = (nx * ny + kBlockThreads - 1) / kBlockThreads;
134  int N = nx * ny;
135  fftShift << <nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, false);
136 
137  cufftHandle plan;
138  cufftResult result;
139  // fft
140  result = cufftPlan2d(&plan, ny, nx, CUFFT_Z2Z);
141  if (result != CUFFT_SUCCESS)
142  {
143  LOG("<FAILED> cufftPlan2d (%d)\n", result);
144  return;
145  };
146 
147  if (direction == -1)
148  result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_FORWARD);
149  else
150  result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_INVERSE);
151 
152  if (result != CUFFT_SUCCESS)
153  {
154  LOG("<FAILED> cufftExecZ2Z (%d)\n", result);
155  return;
156  }
157 
158  if (cudaDeviceSynchronize() != cudaSuccess) {
159  LOG("<FAILED> cudaDeviceSynchronize\n");
160  return;
161  }
162 
163  fftShift << < nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, bNormalized);
164 
165  cufftDestroy(plan);
166  }
167 }
168 
169 #endif // !__ophKernel_cuh__