27194b56
Pavel Govyadinov
major bug fixes, ...
|
1
2
3
4
5
6
7
8
9
|
#ifndef STIM_FILTER_H
#define STIM_FILTER_H
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <stdio.h>
#include <stim/visualization/colormap.h>
#include <sstream>
|
9474c283
Pavel Govyadinov
removed some unus...
|
10
|
#include <stim/math/constants.h>
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
11
12
13
14
15
16
17
|
#include <stim/cuda/cudatools/devices.h>
#include <stim/cuda/cudatools/threads.h>
#include <stim/cuda/cuda_texture.cuh>
#include <stim/cuda/ivote.cuh>
#include <stim/cuda/arraymath.cuh>
#define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
namespace stim
{
namespace cuda
{
float* gpuLoG;
float* LoG;
float* res;
float* centers;
stim::cuda::cuda_texture tx;
void initArray(int DIM_X, int DIM_Y, int kl)
{
LoG = (float*) malloc(kl*kl*sizeof(float));
HANDLE_ERROR(
cudaMalloc( (void**) &gpuLoG, kl*kl*sizeof(float))
);
// checkCUDAerrors("Memory Allocation, LoG");
HANDLE_ERROR(
cudaMalloc( (void**) &res, DIM_Y*DIM_X*sizeof(float))
);
HANDLE_ERROR(
cudaMalloc( (void**) ¢ers, DIM_Y*DIM_X*sizeof(float))
);
// checkCUDAerrors("Memory Allocation, Result");
}
|
c37611a6
Pavel Govyadinov
removed the time ...
|
50
|
void cleanUP()
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
{
HANDLE_ERROR(
cudaFree(gpuLoG)
);
HANDLE_ERROR(
cudaFree(res)
);
HANDLE_ERROR(
cudaFree(centers)
);
free(LoG);
}
void
filterKernel (float kl, float sigma, float *LoG)
{
float t = 0.0;
float kr = kl/2;
float x, y;
int idx;
for(int i = 0; i < kl; i++){
for(int j = 0; j < kl; j++){
idx = j*kl+i;
x = i - kr - 0.5;
y = j - kr - 0.5;
|
9474c283
Pavel Govyadinov
removed some unus...
|
76
|
LoG[idx] = (-1.0/PI/powf(sigma, 4))* (1 - (powf(x,2)+powf(y,2))/2.0/powf(sigma, 2))
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
*expf(-(powf(x,2)+powf(y,2))/2/powf(sigma,2));
t +=LoG[idx];
}
}
for(int i = 0; i < kl*kl; i++)
{
LoG[i] = LoG[i]/t;
}
}
//Shared memory would be better.
__global__
void
applyFilter(cudaTextureObject_t texIn, unsigned int DIM_X, unsigned int DIM_Y, int kr, int kl, float *res, float* gpuLoG){
//R = floor(size/2)
//THIS IS A NAIVE WAY TO DO IT, and there is a better way)
__shared__ float shared[7][7];
int x = blockIdx.x;
int y = blockIdx.y;
int xi = threadIdx.x;
int yi = threadIdx.y;
|
9474c283
Pavel Govyadinov
removed some unus...
|
101
|
// float val = 0;
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
102
103
104
105
106
107
108
109
110
111
112
113
|
float tu = (x-kr+xi)/(float)DIM_X;
float tv = (y-kr+yi)/(float)DIM_Y;
shared[xi][yi] = gpuLoG[yi*kl+xi]*(255.0-(float)tex2D<unsigned char>(texIn, tu, tv));
__syncthreads();
//x = max(0,x);
//x = min(x, width-1);
//y = max(y, 0);
//y = min(y, height - 1);
int idx = y*DIM_X+x;
|
9474c283
Pavel Govyadinov
removed some unus...
|
114
|
// int k_idx;
|
27194b56
Pavel Govyadinov
major bug fixes, ...
|
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
{
__syncthreads();
if (xi < step)
{
shared[xi][yi] += shared[xi + step][yi];
}
__syncthreads();
}
__syncthreads();
for(unsigned int step = blockDim.y/2; step >= 1; step >>= 1)
{
__syncthreads();
if(yi < step)
{
shared[xi][yi] += shared[xi][yi + step];
}
__syncthreads();
}
__syncthreads();
if(xi == 0 && yi == 0)
res[idx] = shared[0][0];
}
extern "C"
float *
get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, float threshold)
{
tx.SetTextureCoordinates(1);
tx.SetAddressMode(1, 3);
tx.MapCudaTexture(texbufferID, texType);
float* result = (float*) malloc(DIM_X*DIM_Y*sizeof(float));
initArray(DIM_X, DIM_Y, sizeK);
filterKernel(sizeK, sigma, LoG);
cudaMemcpy(gpuLoG, LoG, sizeK*sizeK*sizeof(float), cudaMemcpyHostToDevice);
dim3 numBlocks(DIM_X, DIM_Y);
dim3 threadsPerBlock(sizeK, sizeK);
applyFilter <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), DIM_X, DIM_Y, floor(sizeK/2), sizeK, res, gpuLoG);
stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
cudaDeviceSynchronize();
cudaMemcpy(result, centers, DIM_X*DIM_Y*sizeof(float), cudaMemcpyDeviceToHost);
tx.UnmapCudaTexture();
cleanUP();
return result;
}
}
}
#endif
|