Openholo  v4.2
Open Source Digital Holographic Library
ophKernel.cu
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 #include "ophKernel.cuh"
46 
47 __global__ void fftShift(int N, int nx, int ny, cufftDoubleComplex* input, cufftDoubleComplex* output, bool bNormalized)
48 {
49  int tid = threadIdx.x + blockIdx.x * blockDim.x;
50 
51  double normalF = 1.0;
52  if (bNormalized == true)
53  normalF = nx * ny;
54 
55  int half_nx = nx >> 1;
56  int half_ny = ny >> 1;
57  unsigned int size = blockDim.x * gridDim.x;
58 
59  while (tid < N)
60  {
61  int i = tid % nx;
62  int j = tid / nx;
63 
64  int ti = i - half_nx; if (ti < 0) ti += nx;
65  int tj = j - half_ny; if (tj < 0) tj += ny;
66 
67  int oindex = tj * nx + ti;
68 
69 
70  output[tid].x = input[oindex].x / normalF;
71  output[tid].y = input[oindex].y / normalF;
72 
73  tid += size;
74  }
75 }
76 
77 __global__ void fftShiftf(int N, int nx, int ny, cuFloatComplex* input, cuFloatComplex* output, bool bNormalized)
78 {
79  int tid = threadIdx.x + blockIdx.x * blockDim.x;
80 
81  float normalF = 1.0;
82  if (bNormalized == true)
83  normalF = nx * ny;
84 
85  int half_nx = nx >> 1;
86  int half_ny = ny >> 1;
87 
88  unsigned int size = blockDim.x * gridDim.x;
89 
90  while (tid < N)
91  {
92  int i = tid % nx;
93  int j = tid / nx;
94 
95  int ti = i - half_nx; if (ti < 0) ti += nx;
96  int tj = j - half_ny; if (tj < 0) tj += ny;
97 
98  int oindex = tj * nx + ti;
99 
100 
101  output[tid].x = input[oindex].x / normalF;
102  output[tid].y = input[oindex].y / normalF;
103 
104  tid += size;
105  }
106 }
107 
108 __device__ void exponent_complex(cuDoubleComplex* val)
109 {
110  double exp_val = exp(val->x);
111  //double re = val->x;
112  double im = val->y;
113 
114  val->x = exp_val * cos(im);
115  val->y = exp_val * sin(im);
116 }
117 
118 extern "C"
119 {
120  void cudaFFT(CUstream_st* stream, int nx, int ny, cufftDoubleComplex* in_field, cufftDoubleComplex* output_field, int direction, bool bNormalized)
121  {
122  unsigned int nblocks = (nx * ny + kBlockThreads - 1) / kBlockThreads;
123  int N = nx * ny;
124  fftShift << <nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, false);
125 
126  cufftHandle plan;
127  cufftResult result;
128  // fft
129  result = cufftPlan2d(&plan, ny, nx, CUFFT_Z2Z);
130  if (result != CUFFT_SUCCESS)
131  {
132  LOG("<FAILED> cufftPlan2d (%d)\n", result);
133  return;
134  };
135 
136  if (direction == -1)
137  result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_FORWARD);
138  else
139  result = cufftExecZ2Z(plan, output_field, in_field, CUFFT_INVERSE);
140 
141  if (result != CUFFT_SUCCESS)
142  {
143  LOG("<FAILED> cufftExecZ2Z (%d)\n", result);
144  return;
145  }
146 
147  if (cudaDeviceSynchronize() != cudaSuccess) {
148  LOG("<FAILED> cudaDeviceSynchronize\n");
149  return;
150  }
151 
152  fftShift << < nblocks, kBlockThreads, 0, stream >> > (N, nx, ny, in_field, output_field, bNormalized);
153 
154  cufftDestroy(plan);
155  }
156 }