-
Notifications
You must be signed in to change notification settings - Fork 31
/
TVMin3.cc
184 lines (165 loc) · 6.28 KB
/
TVMin3.cc
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
181
182
183
184
#include <TVMin3.h>
#include <vw/Image/Filter.h>
#include <vw/Image/Statistics.h>
#include <vw/FileIO.h>
using namespace vw;
void partial_x( ImageView<float> const& input,
ImageView<float> & output ) {
output.set_size(input.cols(), input.rows());
fill(output,0);
crop(output, BBox2i(0,0,input.cols()-1,input.rows()))
+= crop(input, BBox2i(1,0,input.cols()-1,input.rows()));
crop(output, BBox2i(0,0,input.cols()-1,input.rows()))
-= crop(input, BBox2i(0,0,input.cols()-1,input.rows()));
}
void partial_y( ImageView<float> const& input,
ImageView<float> & output ) {
output.set_size(input.cols(), input.rows());
fill(output,0);
crop(output, BBox2i(0,0,input.cols(),input.rows()-1))
+= crop(input, BBox2i(0,1,input.cols(),input.rows()-1));
crop(output, BBox2i(0,0,input.cols(),input.rows()-1))
-= crop(input, BBox2i(0,0,input.cols(),input.rows()-1));
}
void stereo::divergence( ImageView<float> const& input_x,
ImageView<float> const& input_y,
ImageView<float> & output ) {
output.set_size(input_x.cols(), input_x.rows());
fill(output,0);
crop(output, BBox2i(1,0,input_x.cols()-1,input_x.rows())) -= crop(input_x, BBox2i(0,0,input_x.cols()-1,input_x.rows()));
crop(output, BBox2i(0,0,input_x.cols()-1,input_x.rows())) += crop(input_x, BBox2i(0,0,input_x.cols()-1,input_x.rows()));
crop(output, BBox2i(0,1,input_x.cols(),input_x.rows()-1)) -= crop(input_y, BBox2i(0,0,input_x.cols(),input_x.rows()-1));
crop(output, BBox2i(0,0,input_x.cols(),input_x.rows()-1)) += crop(input_y, BBox2i(0,0,input_x.cols(),input_x.rows()-1));
//output = -diff_x - diff_y;
//output = input_x + input_y;
}
void stereo::gradient( ImageView<float> const& input,
ImageView<float> & output_x,
ImageView<float> & output_y) {
partial_x(input, output_x);
partial_y(input, output_y);
}
float calc_energy_ROF( ImageView<float> const& input,
ImageView<float> const& ref,
float lambda ) {
ImageView<float> dx, dy;
stereo::gradient(input, dx, dy);
float e_reg = sum_of_pixel_values(sqrt(dx * dx + dy * dy));
float e_data = 0.5 * lambda * sum_of_pixel_values((input - ref) * (input - ref));
return e_reg + e_data;
}
void stereo::ROF( ImageView<float> const& input,
float lambda, int iterations,
float sigma, float tau, // These are gradient step sizes
ImageView<float> & output ) {
// Allocate space for p, our hidden variable and u our output.
ImageView<float> p_x(input.cols(), input.rows()),
p_y(input.cols(), input.rows());
output.set_size(input.cols(), input.rows());
output = copy(input);
ImageView<float> grad_u_x, grad_u_y;
ImageView<float> div_p;
stereo::gradient(output, p_x, p_y);
for ( int i = 0; i < iterations; i++ ) {
// Eqn 136
stereo::gradient(output, grad_u_x, grad_u_y);
p_x += sigma * grad_u_x;
p_y += sigma * grad_u_y;
// Eqn 137 (project to a unit sphere
for (int j = 0; j < p_x.rows(); j++ ) {
for (int i = 0; i < p_x.cols(); i++ ) {
float mag =
std::max(1.0, sqrt(p_x(i,j)*p_x(i,j) +
p_y(i,j)*p_y(i,j)));
p_x(i,j) /= mag;
p_y(i,j) /= mag;
}
}
//std::cout << "g x: " << grad_u_x(10,10) << " y: " << grad_u_y(10,10) << " p x: " << p_x(10,10) << " y: " << p_y(10,10) << " " << output(10, 10) << std::endl;
// Eqn 139
stereo::divergence(p_x, p_y, div_p);
output += tau * div_p + tau * lambda * input;
// Eqn 139 (pt2)
output /= (1 + tau * lambda);
//std::cout << "Energy: " << calc_energy_ROF(output, input, lambda) << std::endl;
}
}
void stereo::HuberROF( ImageView<float> const& input,
float lambda, int iterations,
float alpha, // Huber threshold coeff,
float sigma, float tau, // Gradient step sizes
ImageView<float> & output ) {
// Allocate space for p, our hidden variable and u our output.
ImageView<float> p_x(input.cols(), input.rows()),
p_y(input.cols(), input.rows());
output.set_size(input.cols(), input.rows());
output = copy(input);
ImageView<float> grad_u_x, grad_u_y;
ImageView<float> div_p;
gradient(output, p_x, p_y);
for ( int i = 0; i < iterations; i++ ) {
// Eqn 156 & 157
gradient(output, grad_u_x, grad_u_y);
p_x += sigma * grad_u_x;
p_y += sigma * grad_u_y;
p_x /= (1 + sigma * alpha);
p_y /= (1 + sigma * alpha);
// Eqn 158
for (int j = 0; j < p_x.rows(); j++ ) {
for (int i = 0; i < p_x.cols(); i++ ) {
float mag =
std::max(1.0, sqrt(p_x(i,j)*p_x(i,j) +
p_y(i,j)*p_y(i,j)));
p_x(i,j) /= mag;
p_y(i,j) /= mag;
}
}
// Eqn 160
divergence(p_x, p_y, div_p);
output += tau * div_p + tau * lambda * input;
output /= (1 + tau * lambda);
}
}
void stereo::ROF_TVL1( ImageView<float> const& input,
float lambda, int iterations,
float sigma, float tau, // Gradient step sizes
ImageView<float> & output ) {
// Allocate space for p, our hidden variable and u our output.
ImageView<float> p_x(input.cols(), input.rows()),
p_y(input.cols(), input.rows());
output.set_size(input.cols(), input.rows());
output = copy(input);
ImageView<float> grad_u_x, grad_u_y;
ImageView<float> div_p, q;
stereo::gradient(output, p_x, p_y);
q = copy(input);
for ( int i = 0; i < iterations; i++ ) {
// Eqn 185
gradient(output, grad_u_x, grad_u_y);
p_x += sigma * grad_u_x;
p_y += sigma * grad_u_y;
// Eqn 186
for (int j = 0; j < p_x.rows(); j++ ) {
for (int i = 0; i < p_x.cols(); i++ ) {
float mag =
std::max(1.0, sqrt(p_x(i,j)*p_x(i,j) +
p_y(i,j)*p_y(i,j)));
p_x(i,j) /= mag;
p_y(i,j) /= mag;
}
}
// Eqn 188
q += sigma * lambda * (output - input);
// Eqn 189
for (int j = 0; j < q.rows(); j++ ) {
for (int i = 0; i < q.cols(); i++ ) {
float mag =
std::max(1.0, fabs(q(i,j)));
q(i,j) /= mag;
}
}
// Eqn 191
stereo::divergence(p_x, p_y, div_p);
output += tau * div_p - tau * lambda * q;
}
}