PICCANTE  0.4
The hottest HDR imaging library!
mitsunaga_nayar_crf.hpp
Go to the documentation of this file.
1 /*
2 
3 PICCANTE
4 The hottest HDR imaging library!
5 http://vcg.isti.cnr.it/piccante
6 
7 Copyright (C) 2014-2016
8 Visual Computing Laboratory - ISTI CNR
9 http://vcg.isti.cnr.it
10 Library author: Francesco Banterle
11 This file author: Giorgio Marcias
12 
13 This Source Code Form is subject to the terms of the Mozilla Public
14 License, v. 2.0. If a copy of the MPL was not distributed with this
15 file, You can obtain one at http://mozilla.org/MPL/2.0/.
16 
17 */
18 
19 #ifndef PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP
20 #define PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP
21 
22 #include<algorithm>
23 #include<limits>
24 #include<vector>
25 
26 #include "../base.hpp"
27 
28 #ifndef PIC_DISABLE_EIGEN
29 
30 #ifndef PIC_EIGEN_NOT_BUNDLED
31  #include "../externals/Eigen/LU"
32 #else
33  #include <Eigen/LU>
34 #endif
35 
36 #endif
37 
38 namespace pic {
39 
52 PIC_INLINE float MitsunagaNayarClassic(int *samples, const std::size_t nSamples, const std::vector<float> &exposures,
53  std::vector<float> &coefficients, const bool computeRatios, std::vector<float> &R,
54  const float eps, const std::size_t max_iterations)
55 {
56 #ifndef PIC_DISABLE_EIGEN
57  float eval, val, tmp1, tmp2;
58  const std::size_t Q = exposures.size();
59  const std::size_t N = coefficients.size() - 1;
60 
61  const float Mmax = 1.f;
62 
63  for (float &_c : coefficients) {
64  _c = 0.f;
65  }
66 
67  if (!samples || Q < 2 || coefficients.size() < 2) {
68  return std::numeric_limits<float>::infinity();
69  }
70 
71  R.assign(Q < 2 ? 0 : Q - 1, 1.f);
72  for (int q = 0; q < R.size(); ++q) {
73  R[q] = exposures[q] / exposures[q+1];
74  }
75 
76  //Check valid samples
77  std::vector<std::vector<float>> g(nSamples, std::vector<float>(Q, 0.f));
78  std::size_t P = 0;
79  for (std::size_t p = 0; p < nSamples; ++p) {
80  bool valid = false;
81  for (std::size_t q = 0; q < Q-1; ++q) {
82  if (samples[p * Q + q] >= 0 && samples[p * Q + q+1] >= 0) {
83  valid = true;
84  break;
85  }
86  }
87  if (valid) {
88  for (std::size_t q = 0; q < Q; ++q) {
89  if (samples[p * Q + q] >= 0) {
90  g[P][q] = samples[p * Q + q] / 255.f;
91  } else {
92  g[P][q] = -1.f;
93  }
94  }
95  ++P;
96  }
97  }
98  g.resize(P);
99 
100  if (g.empty()) {
101  return std::numeric_limits<float>::infinity();
102  }
103 
104  //Precompute test with exponentials
105  std::vector<Eigen::VectorXf> test;
106  if (computeRatios) {
107  test.assign(256, Eigen::VectorXf::Zero(N+1));
108  for (std::size_t i = 0; i < 256; ++i) {
109  test[i][0] = 1.f;
110  for (std::size_t n = 1; n <= N; ++n) {
111  test[i][n] = (i / 255.f) * test[i][n-1];
112  }
113  }
114  }
115 
116  //Precompute M with exponentials
117  std::vector<std::vector<std::vector<float>>> M(P,
118  std::vector<std::vector<float>>(Q,
119  std::vector<float>(N+1, 0.f)));
120  for (std::size_t p = 0; p < P; ++p) {
121  for (std::size_t q = 0; q < Q; ++q) {
122  M[p][q][0] = 1.f;
123  for (std::size_t n = 1; n <= N; ++n) {
124  M[p][q][n] = g[p][q] * M[p][q][n-1];
125  }
126  }
127  }
128 
129  std::vector<std::vector<std::vector<float>>> d(P,
130  std::vector<std::vector<float>>(Q-1,
131  std::vector<float>(N+1, 1.f)));
132  Eigen::MatrixXf A = Eigen::MatrixXf::Zero(N, N);
133  Eigen::VectorXf x, b = Eigen::VectorXf::Zero(N);
134  Eigen::VectorXf c(N+1), prev_c = Eigen::VectorXf::Zero(N+1);
135 
136  std::size_t iter = 0;
137 
138  do {
139  //Compute d
140  for (std::size_t p = 0; p < P; ++p) {
141  for (std::size_t q = 0; q < Q-1; ++q) {
142  if (g[p][q] >= 0.f && g[p][q+1] >= 0.f) {
143  for (std::size_t n = 0; n <= N; ++n) {
144  d[p][q][n] = M[p][q][n] - R[q] * M[p][q+1][n];
145  }
146  } else {
147  d[p][q].assign(N+1, 0.f);
148  }
149  }
150  }
151 
152  //Build the matrix A of the linear system
153  A.setZero(N, N);
154  for (std::size_t i = 0; i < N; ++i) {
155  for (std::size_t j = 0; j < N; ++j) {
156  for (std::size_t p = 0; p < P; ++p) {
157  for (std::size_t q = 0; q < Q - 1; ++q) {
158  A(i, j) += d[p][q][i] * (d[p][q][j] - d[p][q][N]);
159  }
160  }
161  }
162  }
163 
164  //Build the vector of knowns b
165  b.setZero(N);
166  for (std::size_t i = 0; i < N; ++i) {
167  for (std::size_t p = 0; p < P; ++p) {
168  for (std::size_t q = 0; q < Q - 1; ++q) {
169  b(i) -= Mmax * d[p][q][i] * d[p][q][N];
170  }
171  }
172  }
173 
174  //Solve the linear system
175  x = A.partialPivLu().solve(b);
176  c << x, Mmax - x.sum();
177 
178  if (computeRatios) {
179  //Evaluate approximation increment
180  eval = std::numeric_limits<float>::lowest();
181  for (const Eigen::VectorXf &_M : test) {
182  val = std::abs((c - prev_c).dot(_M));
183  if (val > eval) {
184  eval = val;
185  }
186  }
187 
188  //Update R
189  for (std::size_t q = 0; q < Q-1; ++q) {
190  R[q] = 0.f;
191  tmp1 = 0.f;
192  tmp2 = 0.f;
193  for (std::size_t p = 0; p < P; ++p) {
194  if (g[p][q] >= 0.f && g[p][q+1] >= 0.f) {
195  for (std::size_t n = 0; n <= N; ++n) {
196  tmp1 += c[n] * M[p][q][n];
197  tmp2 += c[n] * M[p][q+1][n];
198  }
199  }
200  }
201  R[q] += tmp1 / tmp2;
202  }
203 
204  ++iter;
205  }
206  } while (computeRatios && eval > eps && iter < max_iterations);
207 
208  for (std::size_t n = 0; n <= N; ++n) {
209  coefficients[n] = c[n];
210  }
211 
212  //Evaluate error
213  eval = 0.f;
214  for (std::size_t q = 0; q < Q-1; ++q) {
215  for (std::size_t p = 0; p < P; ++p) {
216  if (g[p][q] >= 0.f && g[p][q+1] >= 0.f) {
217  val = 0.f;
218  for (std::size_t n = 0; n <= N; ++n) {
219  val += coefficients[n] * (M[p][q][n] - R[q] * M[p][q+1][n]);
220  }
221  eval += val * val;
222  }
223  }
224  }
225 
226  return eval;
227 #else
228  return -1.0f;
229 #endif
230 }
231 
244 PIC_INLINE float MitsunagaNayarFull(int *samples, const std::size_t nSamples, const std::vector<float> &exposures,
245  std::vector<float> &coefficients, bool computeRatios, std::vector<std::vector<float>> &R,
246  const float eps, const std::size_t max_iterations)
247 {
248 #ifndef PIC_DISABLE_EIGEN
249  float eval, val, tmp1, tmp2;
250  const std::size_t Q = exposures.size();
251  const std::size_t N = coefficients.size() - 1;
252 
253  const float Mmax = 1.f;
254 
255  for (float &_c : coefficients) {
256  _c = 0.f;
257  }
258  R.assign(Q < 2 ? 0 : Q, std::vector<float>(Q < 2 ? 0 : Q, 1.f));
259  for (int q1 = 0; q1 < R.size(); ++q1) {
260  for (int q2 = 0; q2 < R[q1].size(); ++q2) {
261  if (q2 == q1) {
262  R[q1][q2] = 1.f;
263  } else {
264  R[q1][q2] = exposures[q1] / exposures[q2];
265  }
266  }
267  }
268  if (!samples || Q < 2 || coefficients.size() < 2) {
269  return std::numeric_limits<float>::infinity();
270  }
271 
272  //Check valid samples
273  std::vector<std::vector<float>> g(nSamples, std::vector<float>(Q, 0.f));
274  std::size_t P = 0;
275  for (std::size_t p = 0; p < nSamples; ++p) {
276  std::size_t valid = 0;
277  for (std::size_t q = 0; q < Q; ++q) {
278  if (samples[p * Q + q] >= 0) {
279  ++valid;
280  }
281  }
282  if (valid > 1) {
283  for (std::size_t q = 0; q < Q; ++q) {
284  if (samples[p * Q + q] >= 0) {
285  g[P][q] = samples[p * Q + q] / 255.f;
286  } else {
287  g[P][q] = -1.f;
288  }
289  }
290  ++P;
291  }
292  }
293  g.resize(P);
294 
295  if (g.empty()) {
296  return std::numeric_limits<float>::infinity();
297  }
298 
299  //Precompute test with exponentials
300  std::vector<Eigen::VectorXf> test;
301  if (computeRatios) {
302  test.assign(256, Eigen::VectorXf::Zero(N+1));
303  for (std::size_t i = 0; i < 256; ++i) {
304  test[i][0] = 1.f;
305  for (std::size_t n = 1; n <= N; ++n) {
306  test[i][n] = (i / 255.f) * test[i][n-1];
307  }
308  }
309  }
310 
311  //Precompute M with exponentials
312  std::vector<std::vector<std::vector<float>>> M(P,
313  std::vector<std::vector<float>>(Q,
314  std::vector<float>(N+1, 0.f)));
315  for (std::size_t p = 0; p < P; ++p) {
316  for (std::size_t q = 0; q < Q; ++q) {
317  M[p][q][0] = 1.f;
318  for (std::size_t n = 1; n <= N; ++n) {
319  M[p][q][n] = g[p][q] * M[p][q][n-1];
320  }
321  }
322  }
323 
324  std::vector<std::vector<std::vector<std::vector<float>>>> d(P,
325  std::vector<std::vector<std::vector<float>>>(Q,
326  std::vector<std::vector<float>>(Q,
327  std::vector<float>(N+1, 1.f))));
328  Eigen::MatrixXf A = Eigen::MatrixXf::Zero(N, N);
329  Eigen::VectorXf x, b = Eigen::VectorXf::Zero(N);
330  Eigen::VectorXf c(N+1), prev_c = Eigen::VectorXf::Zero(N+1);
331 
332  std::size_t iter = 0;
333 
334  do {
335  //Compute d
336  for (std::size_t p = 0; p < P; ++p) {
337  for (std::size_t q1 = 0; q1 < Q; ++q1) {
338  for (std::size_t q2 = 0; q2 < Q; ++q2) {
339  d[p][q1][q2].assign(N+1, 0.f);
340  if (q2 != q1) {
341  for (std::size_t n = 0; n <= N; ++n) {
342  if (g[p][q1] >= 0.f && g[p][q2] >= 0.f) {
343  d[p][q1][q2][n] = M[p][q1][n] - R[q1][q2] * M[p][q2][n];
344  } else {
345  d[p][q1][q2][n] = 0.f;
346  }
347  }
348  }
349  }
350  }
351  }
352 
353  //Build the matrix A of the linear system
354  A.setZero(N, N);
355  for (std::size_t i = 0; i < N; ++i) {
356  for (std::size_t j = 0; j < N; ++j) {
357  for (std::size_t p = 0; p < P; ++p) {
358  for (std::size_t q1 = 0; q1 < Q; ++q1) {
359  for (std::size_t q2 = 0; q2 < Q; ++q2) {
360  if (q2 != q1) {
361  A(i, j) += d[p][q1][q2][i] * (d[p][q1][q2][j] - d[p][q1][q2][N]);
362  }
363  }
364  }
365  }
366  }
367  }
368 
369  //Build the vector of knowns b
370  b.setZero(N);
371  for (std::size_t i = 0; i < N; ++i) {
372  for (std::size_t p = 0; p < P; ++p) {
373  for (std::size_t q1 = 0; q1 < Q; ++q1) {
374  for (std::size_t q2 = 0; q2 < Q; ++q2) {
375  if (q2 != q1) {
376  b(i) -= Mmax * d[p][q1][q2][i] * d[p][q1][q2][N];
377  }
378  }
379  }
380  }
381  }
382 
383  //Solve the linear system
384  x = A.partialPivLu().solve(b);
385  c << x, Mmax - x.sum();
386 
387  if (computeRatios) {
388  //Evaluate approximation increment
389  eval = std::numeric_limits<float>::lowest();
390  for (const Eigen::VectorXf &_M : test) {
391  val = std::abs((c - prev_c).dot(_M));
392  if (val > eval) {
393  eval = val;
394  }
395  }
396 
397  //Update R
398  for (std::size_t q1 = 0; q1 < Q; ++q1) {
399  for (std::size_t q2 = 0; q2 < Q; ++q2) {
400  R[q1][q2] = 0.f;
401  tmp1 = 0.f;
402  tmp2 = 0.f;
403  for (std::size_t p = 0; p < P; ++p) {
404  if (g[p][q1] >= 0.f && g[p][q2] >= 0.f) {
405  for (std::size_t n = 0; n <= N; ++n) {
406  tmp1 += c[n] * M[p][q1][n];
407  tmp2 += c[n] * M[p][q2][n];
408  }
409  }
410  }
411  R[q1][q2] += tmp1 / tmp2;
412  }
413  }
414 
415  ++iter;
416  }
417  } while (computeRatios && eval > eps && iter < max_iterations);
418 
419  for (std::size_t n = 0; n <= N; ++n) {
420  coefficients[n] = c[n];
421  }
422 
423  //Evaluate error
424  eval = 0.f;
425  for (std::size_t q1 = 0; q1 < Q; ++q1) {
426  for (std::size_t q2 = 0; q2 < Q; ++q2) {
427  if (q2 != q1) {
428  for (std::size_t p = 0; p < P; ++p) {
429  if (g[p][q1] >= 0.f && g[p][q2] >= 0.f) {
430  val = 0.f;
431  for (std::size_t n = 0; n <= N; ++n) {
432  val += coefficients[n] * (M[p][q1][n] - R[q1][q2] * M[p][q2][n]);
433  }
434  eval += val * val;
435  }
436  }
437  }
438  }
439  }
440 
441  return eval;
442 #else
443  return -1.0f;
444 #endif
445 }
446 
447 }
448 
449 #endif // PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP
pic::MitsunagaNayarFull
PIC_INLINE float MitsunagaNayarFull(int *samples, const std::size_t nSamples, const std::vector< float > &exposures, std::vector< float > &coefficients, bool computeRatios, std::vector< std::vector< float >> &R, const float eps, const std::size_t max_iterations)
MitsunagaNayarFull computes the inverse CRF of a camera as a polynomial function, using all exposure ...
Definition: mitsunaga_nayar_crf.hpp:244
pic::MitsunagaNayarClassic
PIC_INLINE float MitsunagaNayarClassic(int *samples, const std::size_t nSamples, const std::vector< float > &exposures, std::vector< float > &coefficients, const bool computeRatios, std::vector< float > &R, const float eps, const std::size_t max_iterations)
MitsunagaNayarClassic computes the inverse CRF of a camera as a polynomial function.
Definition: mitsunaga_nayar_crf.hpp:52
PIC_INLINE
#define PIC_INLINE
Definition: base.hpp:33
pic
Definition: bilateral_separation.hpp:25
Generated on Mon Oct 19 2020 08:56:12 for PICCANTE by   doxygen 1.8.14