-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathrt_bf.py
More file actions
117 lines (86 loc) · 4.17 KB
/
rt_bf.py
File metadata and controls
117 lines (86 loc) · 4.17 KB
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
# -*- coding: utf-8 -*-
"""180227_Assign_4 Q3
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/151YjjLC9MhFhp2BQQVlkHpiyi4HofDOy
"""
import numpy as np
import cv2
from scipy import signal, interpolate
from google.colab.patches import cv2_imshow
def flatter_image(img):
flat_image = img.flatten()
edgemin = np.amin(flat_image)
edgemax = np.amax(flat_image)
edgedelta = edgemax - edgemin
return flat_image, edgemin, edgemax,edgedelta
def derived(sigmaspatial, sigmarange, samplespatial, samplerange):
samplespatial = sigmaspatial if (samplespatial is None) else samplespatial
samplerange = sigmarange if (samplerange is None) else samplerange
derivedspatial = sigmaspatial / samplespatial
derivedrange = sigmarange / samplerange
return samplespatial, samplerange, derivedspatial, derivedrange
def sample_dimension(width,height,edgedelta,samplespatial,samplerange,xypadding,zpadding):
sw = int(np.round((width - 1) / samplespatial) + 1 + 2 * xypadding)
sh = int(np.round((height - 1) / samplespatial) + 1 + 2 * xypadding)
sd = int(np.round(edgedelta / samplerange) + 1 + 2 * zpadding)
return sw, sh, sd
def get_dim(image,xgrid,ygrid,edgemin,samplespatial,xypadding,zpadding,samplerange,sw, sh, sd):
dimx = np.around(xgrid / samplespatial) + xypadding
dimy = np.around(ygrid / samplespatial) + xypadding
dimz = np.around((image - edgemin) / samplerange) + zpadding
flatx,flaty,flatz = flat_them(dimx,dimy,dimz)
dim = flatz + flaty * sd + flatx * sw * sd
dim = np.array(dim, dtype=int)
return dimx,dimy,dimz,dim
def flat_them(dimx,dimy,dimz):
flatx = dimx.flatten()
flaty = dimy.flatten()
flatz = dimz.flatten()
return flatx,flaty,flatz
def rt_BF(image, sigmaspatial, sigmarange, samplespatial=None, samplerange=None):
height = image.shape[0]
width = image.shape[1]
samplespatial, samplerange, derivedspatial, derivedrange = derived(sigmaspatial, sigmarange, samplespatial, samplerange)
flat_image, edgemin, edgemax,edgedelta = flatter_image(image)
flat_image = image.flatten()
xypadding = np.round(2 * derivedspatial + 1)
zpadding = np.round(2 * derivedrange + 1)
sw, sh, sd = sample_dimension(width,height,edgedelta,samplespatial,samplerange,xypadding,zpadding)
df = np.zeros(sh * sw * sd)
(ygrid, xgrid) = np.meshgrid(range(width), range(height))
dimx,dimy,dimz,dim= get_dim(image,xgrid,ygrid,edgemin,samplespatial,xypadding,zpadding,samplerange,sw, sh, sd)
df[dim] = flat_image
data = df.reshape(sh, sw, sd)
weights = np.array(data, dtype=bool)
kerneldim = derivedspatial * 2 + 1
kerneldep = 2 * derivedrange * 2 + 1
(gridx, gridy, gridz) = np.meshgrid(range(int(kerneldim)), range(int(kerneldim)), range(int(kerneldep)))
halfkerneldim = np.round(kerneldim / 2)
halfkerneldep = np.round(kerneldep / 2)
gridx -= int(halfkerneldim)
gridy -= int(halfkerneldim)
gridz -= int(halfkerneldep)
gridsqr = ((gridx * gridx + gridy * gridy) / (derivedspatial * derivedspatial)) + ((gridz * gridz) / (derivedrange * derivedrange))
kernel = np.exp(-0.5 * gridsqr)
blurweights = signal.fftconvolve(weights, kernel, mode='same')
blurweights = np.where(blurweights == 0, -2, blurweights)
blurdata = signal.fftconvolve(data, kernel, mode='same')
normalblurdata = blurdata / blurweights
normalblurdata = np.where(blurweights < -1, 0, normalblurdata)
(ygrid, xgrid) = np.meshgrid(range(width), range(height))
dimx = (xgrid / samplespatial) + xypadding
dimy = (ygrid / samplespatial) + xypadding
dimz = (image - edgemin) / samplerange + zpadding
out_img = interpolate.interpn((range(normalblurdata.shape[0]), range(normalblurdata.shape[1]),
range(normalblurdata.shape[2])), normalblurdata, (dimx, dimy, dimz))
return out_img
I = cv2.imread('rome.jpg', cv2.IMREAD_UNCHANGED ).astype(np.float32)/255.0
sigmaspatial = 75
sigmarange = 0.2
out = np.stack([
rt_BF( I[:,:,0], sigmaspatial, sigmarange ),
rt_BF( I[:,:,1], sigmaspatial, sigmarange ),
rt_BF( I[:,:,2], sigmaspatial, sigmarange )], axis=2 )
out = out*255.0
cv2_imshow(out )