19 #ifndef PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP 20 #define PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP 26 #include "../base.hpp" 28 #ifndef PIC_DISABLE_EIGEN 30 #ifndef PIC_EIGEN_NOT_BUNDLED 31 #include "../externals/Eigen/LU" 53 std::vector<float> &coefficients,
const bool computeRatios, std::vector<float> &R,
54 const float eps,
const std::size_t max_iterations)
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;
61 const float Mmax = 1.f;
63 for (
float &_c : coefficients) {
67 if (!samples || Q < 2 || coefficients.size() < 2) {
68 return std::numeric_limits<float>::infinity();
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];
77 std::vector<std::vector<float>> g(nSamples, std::vector<float>(Q, 0.f));
79 for (std::size_t p = 0; p < nSamples; ++p) {
81 for (std::size_t q = 0; q < Q-1; ++q) {
82 if (samples[p * Q + q] >= 0 && samples[p * Q + q+1] >= 0) {
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;
101 return std::numeric_limits<float>::infinity();
105 std::vector<Eigen::VectorXf> test;
107 test.assign(256, Eigen::VectorXf::Zero(N+1));
108 for (std::size_t i = 0; i < 256; ++i) {
110 for (std::size_t n = 1; n <= N; ++n) {
111 test[i][n] = (i / 255.f) * test[i][n-1];
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) {
123 for (std::size_t n = 1; n <= N; ++n) {
124 M[p][q][n] = g[p][q] * M[p][q][n-1];
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);
136 std::size_t iter = 0;
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];
147 d[p][q].assign(N+1, 0.f);
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]);
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];
175 x = A.partialPivLu().solve(b);
176 c << x, Mmax - x.sum();
180 eval = std::numeric_limits<float>::lowest();
181 for (
const Eigen::VectorXf &_M : test) {
182 val = std::abs((c - prev_c).dot(_M));
189 for (std::size_t q = 0; q < Q-1; ++q) {
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];
206 }
while (computeRatios && eval > eps && iter < max_iterations);
208 for (std::size_t n = 0; n <= N; ++n) {
209 coefficients[n] = c[n];
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) {
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]);
245 std::vector<float> &coefficients,
bool computeRatios, std::vector<std::vector<float>> &R,
246 const float eps,
const std::size_t max_iterations)
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;
253 const float Mmax = 1.f;
255 for (
float &_c : coefficients) {
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) {
264 R[q1][q2] = exposures[q1] / exposures[q2];
268 if (!samples || Q < 2 || coefficients.size() < 2) {
269 return std::numeric_limits<float>::infinity();
273 std::vector<std::vector<float>> g(nSamples, std::vector<float>(Q, 0.f));
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) {
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;
296 return std::numeric_limits<float>::infinity();
300 std::vector<Eigen::VectorXf> test;
302 test.assign(256, Eigen::VectorXf::Zero(N+1));
303 for (std::size_t i = 0; i < 256; ++i) {
305 for (std::size_t n = 1; n <= N; ++n) {
306 test[i][n] = (i / 255.f) * test[i][n-1];
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) {
318 for (std::size_t n = 1; n <= N; ++n) {
319 M[p][q][n] = g[p][q] * M[p][q][n-1];
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);
332 std::size_t iter = 0;
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);
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];
345 d[p][q1][q2][n] = 0.f;
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) {
361 A(i, j) += d[p][q1][q2][i] * (d[p][q1][q2][j] - d[p][q1][q2][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) {
376 b(i) -= Mmax * d[p][q1][q2][i] * d[p][q1][q2][N];
384 x = A.partialPivLu().solve(b);
385 c << x, Mmax - x.sum();
389 eval = std::numeric_limits<float>::lowest();
390 for (
const Eigen::VectorXf &_M : test) {
391 val = std::abs((c - prev_c).dot(_M));
398 for (std::size_t q1 = 0; q1 < Q; ++q1) {
399 for (std::size_t q2 = 0; q2 < Q; ++q2) {
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];
411 R[q1][q2] += tmp1 / tmp2;
417 }
while (computeRatios && eval > eps && iter < max_iterations);
419 for (std::size_t n = 0; n <= N; ++n) {
420 coefficients[n] = c[n];
425 for (std::size_t q1 = 0; q1 < Q; ++q1) {
426 for (std::size_t q2 = 0; q2 < Q; ++q2) {
428 for (std::size_t p = 0; p < P; ++p) {
429 if (g[p][q1] >= 0.f && g[p][q2] >= 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]);
449 #endif // PIC_ALGORITHMS_MITSUNAGA_NAYAR_CRF_HPP 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_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
#define PIC_INLINE
Definition: base.hpp:33
Definition: bilateral_separation.hpp:25