MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
patch_sampler.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 "math/defines.h"
11#include "math/matrix.h"
12#include "math/vector.h"
13#include "dmrecon/defines.h"
14#include "dmrecon/mvs_tools.h"
16
18
19PatchSampler::PatchSampler(std::vector<SingleView::Ptr> const& _views,
20 Settings const& _settings,
21 int _x, int _y, float _depth, float _dzI, float _dzJ)
22 : views(_views)
23 , settings(_settings)
24 , midPix(_x,_y)
25 , masterMeanCol(0.f)
26 , depth(_depth)
27 , dzI(_dzI)
28 , dzJ(_dzJ)
29 , success(views.size(), false)
30{
31 SingleView::Ptr refV(views[settings.refViewNr]);
32 mve::ByteImage::ConstPtr masterImg(refV->getScaledImg());
33
34 offset = settings.filterWidth / 2;
35 nrSamples = sqr(settings.filterWidth);
36
37 /* initialize arrays */
38 patchPoints.resize(nrSamples);
39 masterColorSamples.resize(nrSamples);
40 masterViewDirs.resize(nrSamples);
41
42 /* compute patch position and check if it's valid */
44 h[0] = h[1] = offset;
45 topLeft = midPix - h;
46 bottomRight = midPix + h;
47 if (topLeft[0] < 0 || topLeft[1] < 0
48 || bottomRight[0] > masterImg->width()-1
49 || bottomRight[1] > masterImg->height()-1)
50 return;
51
52 /* initialize viewing rays from master view */
53 std::size_t count = 0;
54 for (int j = topLeft[1]; j <= bottomRight[1]; ++j)
55 for (int i = topLeft[0]; i <= bottomRight[0]; ++i)
56 masterViewDirs[count++] = refV->viewRayScaled(i, j);
57
58 /* initialize master color samples and 3d patch points */
59 success[settings.refViewNr] = true;
60 computeMasterSamples();
61 computePatchPoints();
62}
63
64void
65PatchSampler::fastColAndDeriv(std::size_t v, Samples& color, Samples& deriv)
66{
67 success[v] = false;
68 SingleView::Ptr refV = views[settings.refViewNr];
69
70 PixelCoords& imgPos = neighPosSamples[v];
71 imgPos.resize(nrSamples);
72
73 math::Vec3f const& p0 = patchPoints[nrSamples/2];
74 /* compute pixel prints and decide on which MipMap-Level to draw
75 the samples */
76 float mfp = refV->footPrintScaled(p0);
77 float nfp = views[v]->footPrint(p0);
78 if (mfp <= 0.f) {
79 std::cerr << "Error in getFastColAndDerivSamples! "
80 << "footprint in master view: " << mfp << std::endl;
81 throw std::out_of_range("Negative pixel footprint");
82 }
83 if (nfp <= 0.f)
84 return;
85 float ratio = nfp / mfp;
86 int mmLevel = 0;
87 while (ratio < 0.5f) {
88 ++mmLevel;
89 ratio *= 2.f;
90 }
91 mmLevel = views[v]->clampLevel(mmLevel);
92
93 /* compute step size for derivative */
94 math::Vec3f p1(p0 + masterViewDirs[nrSamples/2]);
95 float d = (views[v]->worldToScreen(p1, mmLevel)
96 - views[v]->worldToScreen(patchPoints[12], mmLevel)).norm();
97 if (!(d > 0.f)) {
98 return;
99 }
100 stepSize[v] = 1.f / d;
101
102 /* request according undistorted color image */
103 mve::ByteImage::ConstPtr img(views[v]->getPyramidImg(mmLevel));
104 int const w = img->width();
105 int const h = img->height();
106
107 /* compute image position and gradient direction for each sample
108 point in neighbor image v */
109 std::vector<math::Vec2f> gradDir(nrSamples);
110 for (std::size_t i = 0; i < nrSamples; ++i)
111 {
112 math::Vec3f p0(patchPoints[i]);
113 math::Vec3f p1(patchPoints[i] + masterViewDirs[i] * stepSize[v]);
114 imgPos[i] = views[v]->worldToScreen(p0, mmLevel);
115 // imgPos should be away from image border
116 if (!(imgPos[i][0] > 0 && imgPos[i][0] < w-1 &&
117 imgPos[i][1] > 0 && imgPos[i][1] < h-1)) {
118 return;
119 }
120 gradDir[i] = views[v]->worldToScreen(p1, mmLevel) - imgPos[i];
121 }
122
123 /* draw the samples in the image */
124 color.resize(nrSamples, math::Vec3f(0.f));
125 deriv.resize(nrSamples, math::Vec3f(0.f));
126 colAndExactDeriv(*img, imgPos, gradDir, color, deriv);
127
128 /* normalize the gradient */
129 for (std::size_t i = 0; i < nrSamples; ++i)
130 deriv[i] /= stepSize[v];
131
132 success[v] = true;
133}
134
135float
137{
138 if (neighColorSamples[v].empty())
139 computeNeighColorSamples(v);
140 if (!success[v])
141 return -1.f;
142 assert(success[settings.refViewNr]);
143 math::Vec3f meanY(0.f);
144 for (std::size_t i = 0; i < nrSamples; ++i)
145 meanY += neighColorSamples[v][i];
146 meanY /= (float) nrSamples;
147
148 float sqrDevY = 0.f;
149 float devXY = 0.f;
150 for (std::size_t i = 0; i < nrSamples; ++i)
151 {
152 sqrDevY += (neighColorSamples[v][i] - meanY).square_norm();
153 // Note: master color samples are normalized!
154 devXY += (masterColorSamples[i] - meanX)
155 .dot(neighColorSamples[v][i] - meanY);
156 }
157 float tmp = sqrt(sqrDevX * sqrDevY);
158 assert(!MATH_ISNAN(tmp) && !MATH_ISNAN(devXY));
159 if (tmp > 0)
160 return (devXY / tmp);
161 else
162 return -1.f;
163}
164
165float
166PatchSampler::getNCC(std::size_t u, std::size_t v)
167{
168 if (neighColorSamples[u].empty())
169 computeNeighColorSamples(u);
170 if (neighColorSamples[v].empty())
171 computeNeighColorSamples(v);
172 if (!success[u] || !success[v])
173 return -1.f;
174
175 math::Vec3f meanX(0.f);
176 math::Vec3f meanY(0.f);
177 for (std::size_t i = 0; i < nrSamples; ++i)
178 {
179 meanX += neighColorSamples[u][i];
180 meanY += neighColorSamples[v][i];
181 }
182 meanX /= nrSamples;
183 meanY /= nrSamples;
184
185 float sqrDevX = 0.f;
186 float sqrDevY = 0.f;
187 float devXY = 0.f;
188 for (std::size_t i = 0; i < nrSamples; ++i)
189 {
190 sqrDevX += (neighColorSamples[u][i] - meanX).square_norm();
191 sqrDevY += (neighColorSamples[v][i] - meanY).square_norm();
192 devXY += (neighColorSamples[u][i] - meanX)
193 .dot(neighColorSamples[v][i] - meanY);
194 }
195
196 float tmp = sqrt(sqrDevX * sqrDevY);
197 if (tmp > 0)
198 return (devXY / tmp);
199 else
200 return -1.f;
201}
202
203float
204PatchSampler::getSAD(std::size_t v, math::Vec3f const& cs)
205{
206 if (neighColorSamples[v].empty())
207 computeNeighColorSamples(v);
208 if (!success[v])
209 return -1.f;
210
211 float sum = 0.f;
212 for (std::size_t i = 0; i < nrSamples; ++i) {
213 for (int c = 0; c < 3; ++c) {
214 sum += std::abs(cs[c] * neighColorSamples[v][i][c] -
215 masterColorSamples[i][c]);
216 }
217 }
218 return sum;
219}
220
221float
222PatchSampler::getSSD(std::size_t v, math::Vec3f const& cs)
223{
224 if (neighColorSamples[v].empty())
225 computeNeighColorSamples(v);
226 if (!success[v])
227 return -1.f;
228
229 float sum = 0.f;
230 for (std::size_t i = 0; i < nrSamples; ++i)
231 {
232 for (int c = 0; c < 3; ++c)
233 {
234 float diff = cs[c] * neighColorSamples[v][i][c] -
235 masterColorSamples[i][c];
236 sum += diff * diff;
237 }
238 }
239 return sum;
240}
241
244{
245 std::size_t right = nrSamples/2 + offset;
246 std::size_t left = nrSamples/2 - offset;
247 std::size_t top = offset;
248 std::size_t bottom = nrSamples - 1 - offset;
249
250 math::Vec3f a(patchPoints[right] - patchPoints[left]);
251 math::Vec3f b(patchPoints[top] - patchPoints[bottom]);
252 math::Vec3f normal(a.cross(b));
253 normal.normalize();
254
255 return normal;
256}
257
258void
259PatchSampler::update(float newDepth, float newDzI, float newDzJ)
260{
261 success.clear();
262 success.resize(views.size(), false);
263 depth = newDepth;
264 dzI = newDzI;
265 dzJ = newDzJ;
266 success[settings.refViewNr] = true;
267 computePatchPoints();
268 neighColorSamples.clear();
269 neighDerivSamples.clear();
270 neighPosSamples.clear();
271}
272
273void
274PatchSampler::computePatchPoints()
275{
276 SingleView::Ptr refV = views[settings.refViewNr];
277
278 unsigned int count = 0;
279 for (int j = topLeft[1]; j <= bottomRight[1]; ++j)
280 {
281 for (int i = topLeft[0]; i <= bottomRight[0]; ++i)
282 {
283 float tmpDepth = depth + (i - midPix[0]) * dzI +
284 (j - midPix[1]) * dzJ;
285 if (tmpDepth <= 0.f)
286 {
287 success[settings.refViewNr] = false;
288 return;
289 }
290 patchPoints[count] = refV->camPos + tmpDepth *
291 masterViewDirs[count];
292 ++count;
293 }
294 }
295}
296
297void
298PatchSampler::computeMasterSamples()
299{
300 SingleView::Ptr refV = views[settings.refViewNr];
301 mve::ByteImage::ConstPtr img(refV->getScaledImg());
302
303 /* draw color samples from image and compute mean color */
304 std::size_t count = 0;
305 std::vector<math::Vec2i> imgPos(nrSamples);
306 for (int j = topLeft[1]; j <= bottomRight[1]; ++j)
307 for (int i = topLeft[0]; i <= bottomRight[0]; ++i)
308 {
309 imgPos[count][0] = i;
310 imgPos[count][1] = j;
311 ++count;
312 }
313 getXYZColorAtPix(*img, imgPos, &masterColorSamples);
314
315 masterMeanCol = 0.f;
316 for (std::size_t i = 0; i < nrSamples; ++i)
317 for (int c = 0; c < 3; ++c)
318 {
319 assert(masterColorSamples[i][c] >= 0 &&
320 masterColorSamples[i][c] <= 1);
321 masterMeanCol += masterColorSamples[i][c];
322 }
323
324 masterMeanCol /= 3.f * nrSamples;
325 if (masterMeanCol < 0.01f || masterMeanCol > 0.99f) {
326 success[settings.refViewNr] = false;
327 return;
328 }
329
330 meanX.fill(0.f);
331
332 /* normalize master samples so that average intensity over all
333 channels is 1 and compute mean color afterwards */
334 for (std::size_t i = 0; i < nrSamples; ++i)
335 {
336 masterColorSamples[i] /= masterMeanCol;
337 meanX += masterColorSamples[i];
338 }
339 meanX /= nrSamples;
340 sqrDevX = 0.f;
341
342 /* compute variance (independent from actual mean) */
343 for (std::size_t i = 0; i < nrSamples; ++i)
344 sqrDevX += (masterColorSamples[i] - meanX).square_norm();
345}
346
347void
348PatchSampler::computeNeighColorSamples(std::size_t v)
349{
350 SingleView::Ptr refV = views[settings.refViewNr];
351
352 Samples & color = neighColorSamples[v];
353 PixelCoords & imgPos = neighPosSamples[v];
354 success[v] = false;
355
356 /* compute pixel prints and decide on which MipMap-Level to draw
357 the samples */
358 math::Vec3f const & p0 = patchPoints[nrSamples/2];
359 float mfp = refV->footPrintScaled(p0);
360 float nfp = views[v]->footPrint(p0);
361 if (mfp <= 0.f) {
362 std::cerr << "Error in computeNeighColorSamples! "
363 << "footprint in master view: " << mfp << std::endl;
364 throw std::out_of_range("Negative pixel print");
365 }
366 if (nfp <= 0.f)
367 return;
368 float ratio = nfp / mfp;
369
370 int mmLevel = 0;
371 while (ratio < 0.5f) {
372 ++mmLevel;
373 ratio *= 2.f;
374 }
375 mmLevel = views[v]->clampLevel(mmLevel);
376 mve::ByteImage::ConstPtr img(views[v]->getPyramidImg(mmLevel));
377 int const w = img->width();
378 int const h = img->height();
379
380 color.resize(nrSamples);
381 imgPos.resize(nrSamples);
382
383 for (std::size_t i = 0; i < nrSamples; ++i) {
384 imgPos[i] = views[v]->worldToScreen(patchPoints[i], mmLevel);
385 // imgPos should be away from image border
386 if (!(imgPos[i][0] > 0 && imgPos[i][0] < w-1 &&
387 imgPos[i][1] > 0 && imgPos[i][1] < h-1)) {
388 return;
389 }
390 }
391 getXYZColorAtPos(*img, imgPos, &color);
392 success[v] = true;
393}
394
395
Vector class for arbitrary dimensions and types.
Definition vector.h:87
Vector< T, N > cross(Vector< T, N > const &other) const
Cross product between this and another vector.
Definition vector.h:549
Vector< T, N > & fill(T const &value)
Fills all vector elements with the given value.
Definition vector.h:381
Vector< T, N > & normalize(void)
Normalizes self and returns reference to self.
Definition vector.h:448
std::shared_ptr< Image< T > const > ConstPtr
Definition image.h:43
std::vector< bool > success
math::Vec3f getPatchNormal() const
void fastColAndDeriv(std::size_t v, Samples &color, Samples &deriv)
Draw color samples and derivatives in neighbor view v.
float getSSD(std::size_t v, math::Vec3f const &cs)
Computes the sum of squared differences between reference view and neighbor v with respect to color s...
void update(float newDepth, float newDzI, float newDzJ)
float getFastNCC(std::size_t v)
Compute NCC between reference view and a neighbor view.
float getNCC(std::size_t u, std::size_t v)
Compute NCC between two neighboring views.
float getSAD(std::size_t v, math::Vec3f const &cs)
Computes the sum of absolute differences between reference view and neighbor v with respect to color ...
std::shared_ptr< SingleView > Ptr
Definition single_view.h:31
#define MVS_NAMESPACE_BEGIN
Definition defines.h:18
#define MVS_NAMESPACE_END
Definition defines.h:19
#define MATH_ISNAN(x)
Definition defines.h:104
void colAndExactDeriv(mve::ByteImage const &img, PixelCoords const &imgPos, PixelCoords const &gradDir, Samples &color, Samples &deriv)
interpolate color and derivative at given sample positions
Definition mvs_tools.cc:98
void getXYZColorAtPos(mve::ByteImage const &img, PixelCoords const &imgPos, Samples *color)
interpolate only color at given sample positions
Definition mvs_tools.cc:169
std::vector< math::Vec3f > Samples
Definition defines.h:25
const T sqr(const T &a)
Definition defines.h:31
void getXYZColorAtPix(mve::ByteImage const &img, std::vector< math::Vec2i > const &imgPos, Samples *color)
get color at given pixel positions (no interpolation)
Definition mvs_tools.cc:150
std::vector< math::Vec2f > PixelCoords
Definition defines.h:26
unsigned int filterWidth
Size of the patch is width * width, defaults to 5x5.
Definition settings.h:31
std::size_t refViewNr
The reference view ID to reconstruct.
Definition settings.h:25