MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
patch_optimization.cc
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Ronny Klowsky, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#include "util/strings.h"
11#include "math/algo.h"
12#include "math/defines.h"
13#include "math/matrix.h"
14#include "math/matrix_tools.h"
15#include "math/vector.h"
17#include "dmrecon/settings.h"
18
20
21PatchOptimization::PatchOptimization(
22 std::vector<SingleView::Ptr> const& _views,
23 Settings const& _settings,
24 int _x, // Pixel position
25 int _y,
26 float _depth,
27 float _dzI,
28 float _dzJ,
29 IndexSet const & _globalViewIDs,
30 IndexSet const & _localViewIDs)
31 :
32 views(_views),
33 settings(_settings),
34 midx(_x),
35 midy(_y),
36 depth(_depth),
37 dzI(_dzI),
38 dzJ(_dzJ),
39 sampler(PatchSampler::create(views, settings, midx, midy, depth, dzI, dzJ)),
40 ii(sqr(settings.filterWidth)),
41 jj(sqr(settings.filterWidth)),
42 localVS(views, settings, _globalViewIDs, _localViewIDs, sampler)
43{
44 status.iterationCount = 0;
45 status.optiSuccess = true;
46 status.converged = false;
47
48 if (!sampler->success[settings.refViewNr]) {
49 // Sampler could not be initialized properly
50 status.optiSuccess = false;
51 return;
52 }
53
54 std::size_t count = 0;
55
56 int halfFW = (int) settings.filterWidth / 2;
57 pixel_weight.resize(sampler->getNrSamples());
58 for (int j = -halfFW; j <= halfFW; ++j)
59 for (int i = -halfFW; i <= halfFW; ++i) {
60 ii[count] = i;
61 jj[count] = j;
62 pixel_weight[count] = 1.f;
63 ++count;
64 }
65
66 localVS.performVS();
67 if (!localVS.success) {
68 status.optiSuccess = false;
69 return;
70 }
71
72 // brute force initialize all colorScale entries
73 float masterMeanCol = sampler->getMasterMeanColor();
74 for (std::size_t idx = 0; idx < views.size(); ++idx) {
75 colorScale[idx] = math::Vec3f(1.f / masterMeanCol);
76 }
78}
79
80void
82{
83 if (!settings.useColorScale)
84 return;
85 Samples const & mCol = sampler->getMasterColorSamples();
86 IndexSet const & neighIDs = localVS.getSelectedIDs();
87 IndexSet::const_iterator id;
88 for (id = neighIDs.begin(); id != neighIDs.end(); ++id)
89 {
90 // just copied from old mvs:
91 Samples const & nCol = sampler->getNeighColorSamples(*id);
92 if (!sampler->success[*id])
93 return;
94 for (int c = 0; c < 3; ++c) {
95 // for each color channel
96 float ab = 0.f;
97 float aa = 0.f;
98 for (std::size_t i = 0; i < mCol.size(); ++i) {
99 ab += (mCol[i][c] - nCol[i][c] * colorScale[*id][c]) * nCol[i][c];
100 aa += sqr(nCol[i][c]);
101 }
102 if (std::abs(aa) > 1e-6) {
103 colorScale[*id][c] += ab / aa;
104 if (colorScale[*id][c] > 1e3)
105 status.optiSuccess = false;
106 }
107 else
108 status.optiSuccess = false;
109 }
110 }
111}
112
113float
115{
116 SingleView::Ptr refV = views[settings.refViewNr];
117 if (!status.converged)
118 return 0.f;
119
120 /* Compute mean NCC between reference view and local neighbors,
121 where each NCC has to be higher than acceptance NCC */
122 float meanNCC = 0.f;
123 IndexSet const & neighIDs = localVS.getSelectedIDs();
124 IndexSet::const_iterator id;
125 for (id = neighIDs.begin(); id != neighIDs.end(); ++id) {
126 meanNCC += sampler->getFastNCC(*id);
127 }
128 meanNCC /= neighIDs.size();
129
130 float score = (meanNCC - settings.acceptNCC) /
131 (1.f - settings.acceptNCC);
132
133 /* Compute angle between estimated surface normal and view direction
134 and weight current score with dot product */
135 math::Vec3f viewDir(refV->viewRayScaled(midx, midy));
136 math::Vec3f normal(sampler->getPatchNormal());
137 float dotP = - normal.dot(viewDir);
138 if (dotP < 0.2f) {
139 return 0.f;
140 }
141 return score;
142}
143
144float
146{
147 IndexSet const & neighIDs = localVS.getSelectedIDs();
148 IndexSet::const_iterator id;
149 std::size_t nrSamples = sampler->getNrSamples();
150
151 float norm(0);
152 for (id = neighIDs.begin(); id != neighIDs.end(); ++id)
153 {
154 Samples nCol, nDeriv;
155 sampler->fastColAndDeriv(*id, nCol, nDeriv);
156 if (!sampler->success[*id]) {
157 status.optiSuccess = false;
158 return -1.f;
159 }
160
161 math::Vec3f cs(colorScale[*id]);
162 for (std::size_t i = 0; i < nrSamples; ++i) {
163 norm += pixel_weight[i] * (cs.cw_mult(nDeriv[i])).square_norm();
164 }
165 }
166 return norm;
167}
168
169void
171{
172 if (!localVS.success || !status.optiSuccess) {
173 return;
174 }
175
176 // first four iterations only refine depth
177 while ((status.iterationCount < 4) && (status.optiSuccess)) {
179 ++status.iterationCount;
180 }
181
182 bool viewRemoved = false;
183 bool converged = false;
184 while (status.iterationCount < settings.maxIterations &&
185 localVS.success && status.optiSuccess)
186 {
187 IndexSet const & neighIDs = localVS.getSelectedIDs();
188 IndexSet::const_iterator id;
189
190 std::vector<float> oldNCC;
191 for (id = neighIDs.begin(); id != neighIDs.end(); ++id) {
192 oldNCC.push_back(sampler->getFastNCC(*id));
193 }
194
195 status.optiSuccess = false;
196 if (status.iterationCount % 5 == 4 || viewRemoved) {
199 viewRemoved = false;
200 }
201 else {
203 }
204
205 if (!status.optiSuccess)
206 return;
207 converged = true;
208 std::size_t count = 0;
209 IndexSet toBeReplaced;
210
211 for (id = neighIDs.begin(); id != neighIDs.end(); ++id, ++count)
212 {
213 float ncc = sampler->getFastNCC(*id);
214 if (std::abs(ncc - oldNCC[count]) > settings.minRefineDiff)
215 converged = false;
216 if ((ncc < settings.acceptNCC) ||
217 (status.iterationCount == 14 &&
218 std::abs(ncc - oldNCC[count]) > settings.minRefineDiff))
219 {
220 toBeReplaced.insert(*id);
221 viewRemoved = true;
222 }
223 }
224
225 if (viewRemoved) {
226 localVS.replaceViews(toBeReplaced);
227 if (!localVS.success) {
228 return;
229 }
231 }
232 else if (!status.optiSuccess) {
233 // all views valid but optimization failed -> give up
234 return;
235 }
236 else if (converged) {
237 status.converged = true;
238 return;
239 }
240 ++status.iterationCount;
241 }
242}
243
244float
246{
247 Samples const & mCol = sampler->getMasterColorSamples();
248 std::size_t nrSamples = sampler->getNrSamples();
249 IndexSet const & neighIDs = localVS.getSelectedIDs();
250 IndexSet::const_iterator id;
251 float obj = 0.f;
252 for (id = neighIDs.begin(); id != neighIDs.end(); ++id) {
253 Samples const & nCol = sampler->getNeighColorSamples(*id);
254 if (!sampler->success[*id])
255 return -1.f;
256 math::Vec3f cs(colorScale[*id]);
257 for (std::size_t i = 0; i < nrSamples; ++i) {
258 obj += pixel_weight[i] * (mCol[i] - cs.cw_mult(nCol[i])).square_norm();
259 }
260 }
261 return obj;
262}
263
264void
266{
267 float numerator = 0.f;
268 float denom = 0.f;
269 Samples const & mCol = sampler->getMasterColorSamples();
270 IndexSet const & neighIDs = localVS.getSelectedIDs();
271 IndexSet::const_iterator id;
272 std::size_t nrSamples = sampler->getNrSamples();
273
274 for (id = neighIDs.begin(); id != neighIDs.end(); ++id)
275 {
276 Samples nCol, nDeriv;
277 sampler->fastColAndDeriv(*id, nCol, nDeriv);
278 if (!sampler->success[*id]) {
279 status.optiSuccess = false;
280 return;
281 }
282
283 math::Vec3f cs(colorScale[*id]);
284 for (std::size_t i = 0; i < nrSamples; ++i) {
285 numerator += pixel_weight[i] * (cs.cw_mult(nDeriv[i])).dot
286 (mCol[i] - cs.cw_mult(nCol[i]));
287 denom += pixel_weight[i] * (cs.cw_mult(nDeriv[i])).square_norm();
288 }
289 }
290
291 if (denom > 0) {
292 depth += numerator / denom;
293 sampler->update(depth, dzI, dzJ);
294 if (sampler->success[settings.refViewNr])
295 status.optiSuccess = true;
296 else
297 status.optiSuccess = false;
298 }
299}
300
301void
303{
304 if (!localVS.success) {
305 return;
306 }
307 IndexSet const & neighIDs = localVS.getSelectedIDs();
308 std::size_t nrSamples = sampler->getNrSamples();
309
310 // Solve linear system A*x = b using Moore-Penrose pseudoinverse
311 // Fill matrix ATA and vector ATb:
312 math::Matrix3d ATA(0.f);
313 math::Vec3d ATb(0.f);
314 Samples const & mCol = sampler->getMasterColorSamples();
315 IndexSet::const_iterator id;
316 std::size_t row = 0;
317 for (id = neighIDs.begin(); id != neighIDs.end(); ++id)
318 {
319 Samples nCol, nDeriv;
320 sampler->fastColAndDeriv(*id, nCol, nDeriv);
321 if (!sampler->success[*id]) {
322 status.optiSuccess = false;
323 return;
324 }
325 math::Vec3f cs(colorScale[*id]);
326 for (std::size_t i = 0; i < nrSamples; ++i) {
327 for (int c = 0; c < 3; ++c) {
328 math::Vec3f a_i;
329 a_i[0] = pixel_weight[i] * cs[c] * nDeriv[i][c];
330 a_i[1] = pixel_weight[i] * ii[i] * cs[c] * nDeriv[i][c];
331 a_i[2] = pixel_weight[i] * jj[i] * cs[c] * nDeriv[i][c];
332 float b_i = pixel_weight[i] * (mCol[i][c] - cs[c] * nCol[i][c]);
333 assert(!MATH_ISINF(a_i[0]));
334 assert(!MATH_ISINF(a_i[1]));
335 assert(!MATH_ISINF(a_i[2]));
336 assert(!MATH_ISINF(b_i));
337 ATA(0,0) += a_i[0] * a_i[0];
338 ATA(0,1) += a_i[0] * a_i[1];
339 ATA(0,2) += a_i[0] * a_i[2];
340 ATA(1,1) += a_i[1] * a_i[1];
341 ATA(1,2) += a_i[1] * a_i[2];
342 ATA(2,2) += a_i[2] * a_i[2];
343 ATb += a_i * b_i;
344 ++row;
345 }
346 }
347 }
348 ATA(1,0) = ATA(0,1); ATA(2,0) = ATA(0,2); ATA(2,1) = ATA(1,2);
349 double detATA = math::matrix_determinant(ATA);
350 if (detATA == 0.f) {
351 status.optiSuccess = false;
352 return;
353 }
354 math::Matrix3d ATAinv = math::matrix_inverse(ATA, detATA);
355 math::Vec3f X = ATAinv * ATb;
356
357 // update depth and normal
358 dzI += X[1];
359 dzJ += X[2];
360 depth += X[0];
361 sampler->update(depth, dzI, dzJ);
362 if (sampler->success[settings.refViewNr])
363 status.optiSuccess = true;
364 else
365 status.optiSuccess = false;
366}
367
368
Vector class for arbitrary dimensions and types.
Definition vector.h:87
Vector< T, N > cw_mult(Vector< T, N > const &other) const
Component-wise multiplication with another vector.
Definition vector.h:556
T dot(Vector< T, N > const &other) const
Dot (or scalar) product between self and another vector.
Definition vector.h:542
void replaceViews(IndexSet const &toBeReplaced)
std::shared_ptr< SingleView > Ptr
Definition single_view.h:31
IndexSet const & getSelectedIDs() const
#define MVS_NAMESPACE_BEGIN
Definition defines.h:18
#define MVS_NAMESPACE_END
Definition defines.h:19
#define MATH_ISINF(x)
Definition defines.h:105
Matrix< T, N, N > matrix_inverse(Matrix< T, N, N > const &mat)
Calculates the inverse of the given matrix.
T matrix_determinant(Matrix< T, N, N > const &mat)
Calculates the determinant of the given matrix.
Vector< float, 3 > Vec3f
Definition vector.h:31
std::set< std::size_t > IndexSet
Definition defines.h:24
std::vector< math::Vec3f > Samples
Definition defines.h:25
const T sqr(const T &a)
Definition defines.h:31
float minRefineDiff
Definition settings.h:35
float acceptNCC
Definition settings.h:34
unsigned int filterWidth
Size of the patch is width * width, defaults to 5x5.
Definition settings.h:31
bool useColorScale
Definition settings.h:40
unsigned int maxIterations
Definition settings.h:36
std::size_t refViewNr
The reference view ID to reconstruct.
Definition settings.h:25
std::size_t iterationCount