-
Notifications
You must be signed in to change notification settings - Fork 6
/
flowdepth.py
123 lines (88 loc) · 4.26 KB
/
flowdepth.py
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
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 3 16:29:58 2015
@author: gottfried
"""
import os, sys
import numpy as np
import matplotlib.pyplot as plt
from PyQt4 import QtGui
import flowtools as ft
import warping as warp
import solver as pd
def run_ctf(params, d0, data, imageWidget=None):
levels = (np.floor( (np.log(params['minSize'])-np.log(np.minimum(*d0.shape))) / np.log(params['scalefactor']) ) + 1).astype('int')
levels = np.minimum(levels, params['levels'])
warps = params['warps']
iterations = params['iterations']
L = params['Lambda'] # compute lambda sequence for the levels
Lambda = np.zeros(levels)
for idx,l in enumerate(Lambda):
Lambda[idx] = L*(1.0/params['scalefactor'])**idx
params['Lambda'] = Lambda
dim_dual = 3 # dimension of dual variable
for level in range(levels-1,params['stop_level']-1,-1):
level_sz = np.round(np.array(d0.shape) * params['scalefactor']**level)
factor = np.hstack(( level_sz / np.array(d0.shape), 1 ))
K = data['K'] * (np.ones((3,3))*factor).T; K[2,2] = 1 # scale K matrix according to image size
# individual scaling for x & y
print '--- level %d, size %.1f %.1f' % (level, level_sz[1], level_sz[0])
if level == levels-1:
d = ft.imresize(d0, level_sz) # initialization at coarsest level
p = np.zeros(d.shape + (dim_dual,))
else:
print 'prolongate'
d = ft.imresize(d, level_sz) # prolongate to finer level
ptmp = p.copy()
p = np.zeros(d.shape+(dim_dual,))
for i in range(0, dim_dual):
p[:,:,i] = ft.imresize(ptmp[:,:,i], level_sz)
Xn = ft.backproject(d.shape, K)
L_normal = ft.make_linearOperator(d.shape, Xn, K)
img_scaled = [] # scale all images
for img in data['images']:
img_scaled.append(ft.imresize(img, level_sz))
Iref = img_scaled[0]
Gref = data['G'][0]
for k in range(1,warps+1):
print 'warp', k, 'maxZ=', params['maxz']
warpdata = [] # warp all images
for img,g in zip(img_scaled[1:], data['G'][1:]):
Grel = ft.relative_transformation(Gref, g)
warpdata.append(warp.warp_zeta(Iref, img, Grel, K, ft.to_zeta(d)))
d,p,energy = pd.solve_area_pd(warpdata, d, p, iterations, params, params['Lambda'][level], L_normal, imageWidget)
return d
if __name__ == '__main__':
if QtGui.QApplication.instance()==None:
app=QtGui.QApplication(sys.argv)
plt.close('all')
if 'imageWidget' in locals() and imageWidget is not None:
imageWidget.close()
imageWidget = None
npz = np.load('data.npz')
images = [npz['I1'], npz['I2']]; G = [npz['G1'], npz['G2']]
data = dict(images=images, G=G, K=npz['K'])
d_init = 1.8
I1 = data['images'][0]; I2 = data['images'][1]; K = data['K']
d0 = np.ones(I1.shape)*d_init
minalpha = 0.015 # 3d points must be seen at least under this angle
maxz = ft.camera_baseline(data['G'][0], data['G'][1])/np.arctan(minalpha/2) # compute depth value
params = dict(warps=15,
iterations=30,
scalefactor=0.75,
levels=100,
minSize=48,
stop_level=0,
check=10,
Lambda=0.004,
ref=0,
epsilon = 0.001, # huber epsilon
minz=0.5, # lower threshold for scene depth
maxz=maxz) # maxz: clamp depth to be smaller than this value
imageWidget = ft.MplWidget()
imageWidget.show()
plt.figure('I1'); plt.imshow(I1, cmap='gray'); plt.title('Image 1')
plt.figure('I2'); plt.imshow(I2, cmap='gray'); plt.title('Image 2')
d = run_ctf(params, d0, data, imageWidget)
#plt.figure(); plt.imshow(d); plt.colorbar(); plt.title('result')
plt.show()