-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbilateralFilter.pyx
180 lines (136 loc) · 6.02 KB
/
bilateralFilter.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
50
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
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
174
175
176
177
178
179
180
import cython
from cython.parallel import prange
from cython.parallel import parallel
from libc.math cimport exp
from libc.math cimport sqrt
from math import pi as pi
import numpy as np
ctypedef unsigned char uint8_t
@cython.boundscheck(False) # do not check bounds (faster)
@cython.wraparound(False)
def bilateralFilterFast(uint8_t[:, :, :] img, uint8_t sigma_s, uint8_t sigma_b):
# img shape
cdef int H = img.shape[0]
cdef int W = img.shape[1]
cdef int C = img.shape[2]
# pixel values
cdef uint8_t img_r, img_g, img_b, center_r, center_g, center_b
# img after bilateral filter
cdef float[:, :, :] img_filtered = np.zeros([H, W, 3], dtype=np.float32)
# for sigma_s, optimal size is 2 * pi * sigma_s
cdef int dim = int(2 * pi * sigma_s)
# aim dim to be odd
if dim % 2 == 0:
dim += 1
# dim = 2 * k + 1 -> k = (dim - 1) / 2
cdef int k = int((dim - 1) / 2)
# compute gaussian offline
cdef float[:, :] gaussian = np.zeros((dim, dim), dtype = np.float32)
# norm factor for spatial gaussian
cdef float ws = 1 / (2 * pi * sigma_s ** 2)
# norm factor for brightness gaussian
cdef float wb = 1 / (sqrt(2 * pi) * sigma_b)
cdef float wsb, wt
# iterators
cdef int i, j, h, w, x, y
# compute gaussian (m^2 + n^2) / (2 * sigma_s^2)
for i in range(-k, k+1):
for j in range(-k, k+1):
gaussian[i+k, j+k] = (sqr(i) + sqr(j)) / (2 * sigma_s ** 2)
# openmp for faster computation
with nogil, parallel():
for h in prange(H, schedule = 'guided'):
for w in range(W):
# center pixel in the KxK gaussian kernel
center_r = img[h, w, 0]
center_g = img[h, w, 1]
center_b = img[h, w, 2]
wsb = 0
for i in range(-k, k+1):
for j in range(-k, k+1):
# if out of bounds, it simulates padding the borders with the elements at the edge
x = min(max(0, h + i), H - 1)
y = min(max(0, w + j), W - 1)
img_r = img[x, y, 0]
img_g = img[x, y, 1]
img_b = img[x, y, 2]
# ws * wb * exp(-1/2 ( (m^2+n^2)/sigma_s^2 + ( dif_r^2 + dif_g^2 + dif_b^2)/sigma_b^2))
# (m^2 + n^2)/(2 * sigma_s^2) is the gaussian computed above
wt = ws * wb * exp(- (gaussian[i+k, j+k] + (sqr(img_r - center_r) + sqr(img_g - center_g) + sqr(img_b - center_b)) / (2 * sigma_b**2)))
img_filtered[h, w, 0] += img_r * wt
img_filtered[h, w, 1] += img_g * wt
img_filtered[h, w, 2] += img_b * wt
# add to factor Wsb
wsb += wt
# 1e-9 add for numerical stability
img_filtered[h, w, 0] /= (wsb + 1e-9)
img_filtered[h, w, 1] /= (wsb + 1e-9)
img_filtered[h, w, 2] /= (wsb + 1e-9)
return img_filtered
cdef inline float sqr(float x) nogil:
return x * x
"""
def computeGaussian(int sigma_s):
cdef int dim = int(2 * pi * sigma_s)
cdef int k, i, j
if dim % 2 == 0:
dim += 1
k = int((dim - 1) / 2)
cdef float[:, :] gaussian = np.zeros((dim, dim), dtype = np.float32)
cdef float norm_factor = 1 / (2 * pi * sigma_s ** 2)
for i in range(-k, k+1, 1):
for j in range(-k, k+1, 1):
gaussian[i+k, j+k] = norm_factor * exp(-1/2 * (i**2 + j**2) / (sigma_s ** 2))
return np.asarray(gaussian), k
"""
"""
def bilateralFilterFast(float[:, :, :] img, int sigma_s, int sigma_b):
cdef int H = img.shape[0]
cdef int W = img.shape[1]
cdef int C = img.shape[2]
cdef float[:, :, :] img_filtered = np.zeros([H, W, C], dtype=np.float32)
cdef float center_r, center_g, center_b
cdef float wsb_r, wsb_g, wsb_b
cdef float sum_r, sum_g, sum_b
cdef float value_r, value_g, value_b
cdef float dif_r, dif_g, dif_b
cdef float spatial, tonal_r, tonal_g, tonal_b
cdef float norm_factor = 1 / (2 * pi * sigma_s ** 2)
cdef int h, w, i, j, x, y
cdef int k = int(sqrt(sigma_s) * 3)
# The rest of the code is similar, only now we have to explicitly assign the r, g, and b channels
for h in range(H):
for w in range(W):
sum_r = 0
sum_g = 0
sum_b = 0
wsb = 0
for i in range(-k, k):
for j in range(-k, k):
x = np.clip(h + i, 0, H-1)
y = np.clip(w + j, 0, W-1)
value_r = img[x, y, 0]
value_g = img[x, y, 1]
value_b = img[x, y, 2]
dif_r = center_r - value_r
dif_g = center_g - value_g
dif_b = center_b - value_b
spatial = norm_factor * exp(-1/2 * (i**2 + j**2) / (sigma_s ** 2))
# tonal_r = 1 / (sqrt(2 * pi * sigma_b)) * exp(- 1 / 2 * (dif_r/sigma_b**2))
# tonal_g = 1 / (sqrt(2 * pi * sigma_b)) * exp(- 1 / 2 * (dif_g/sigma_b**2))
# tonal_b = 1 / (sqrt(2 * pi * sigma_b)) * exp(- 1 / 2 * (dif_b/sigma_b**2))
# sum_r += value_r * spatial * tonal_r
# sum_g += value_g * spatial * tonal_g
# sum_b += value_b * spatial * tonal_b
# wsb_r += spatial * tonal_r
# wsb_g += spatial * tonal_g
# wsb_b += spatial * tonal_b
sum_r += spatial * value_r
sum_g += spatial * value_g
sum_b += spatial * value_b
wsb += spatial
img_filtered[h, w, 0] = sum_r / (wsb + 1e-6)
img_filtered[h, w, 1] = sum_g / (wsb + 1e-6)
img_filtered[h, w, 2] = sum_b / (wsb + 1e-6)
return img_filtered
"""