Line data Source code
1 : // Copyright (c) 2018-2022 Made to Order Software Corp. All Rights Reserved
2 : //
3 : // https://snapwebsites.org/project/snapdev
4 : // contact@m2osw.com
5 : //
6 : // This program is free software; you can redistribute it and/or modify
7 : // it under the terms of the GNU General Public License as published by
8 : // the Free Software Foundation; either version 2 of the License, or
9 : // (at your option) any later version.
10 : //
11 : // This program is distributed in the hope that it will be useful,
12 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : // GNU General Public License for more details.
15 : //
16 : // You should have received a copy of the GNU General Public License along
17 : // with this program; if not, write to the Free Software Foundation, Inc.,
18 : // 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 : #pragma once
20 :
21 : // C++ lib
22 : //
23 : #include <cctype>
24 : #include <cmath>
25 : #include <iostream>
26 : #include <sstream>
27 : #include <stdexcept>
28 : #include <vector>
29 :
30 :
31 :
32 : namespace snapdev
33 : {
34 :
35 : // Matrix additions, subtractions, and multiplications can be verified
36 : // using this web page:
37 : // http://www.calcul.com/show/calculator/matrix-multiplication
38 :
39 :
40 :
41 : template<typename T, typename SIZE = std::size_t>
42 143362 : class matrix
43 : {
44 : public:
45 : typedef T value_type;
46 : typedef SIZE size_type;
47 :
48 : // The color weights are used to convert RGB to Luma.
49 : //
50 : // With a factor, it's possible to change the color toward or
51 : // away from the perfect Luma if the color is not already
52 : // a gray color (see the saturation() function.)
53 : //
54 : // Note that these are often referenced to as Luminance Weights.
55 : // The Luminance is what you get on your monitor itself, not with
56 : // linear RGB as we manage in software.
57 : //
58 : // The AVERAGE luma factors are present because they may be useful
59 : // in some situation, but they are definitely wrong and should
60 : // very rarely be used.
61 : //
62 : // see https://en.wikipedia.org/wiki/Luma_%28video%29
63 : // see https://www.opengl.org/archives/resources/code/samples/advanced/advanced97/notes/node140.html
64 : //
65 : static value_type constexpr HDTV_LUMA_RED = 0.2126;
66 : static value_type constexpr HDTV_LUMA_GREEN = 0.7152;
67 : static value_type constexpr HDTV_LUMA_BLUE = 0.0722;
68 :
69 : static value_type constexpr LED_LUMA_RED = 0.212;
70 : static value_type constexpr LED_LUMA_GREEN = 0.701;
71 : static value_type constexpr LED_LUMA_BLUE = 0.087;
72 :
73 : static value_type constexpr CRT_LUMA_RED = 0.3086;
74 : static value_type constexpr CRT_LUMA_GREEN = 0.6094;
75 : static value_type constexpr CRT_LUMA_BLUE = 0.0820;
76 :
77 : static value_type constexpr NTSC_LUMA_RED = 0.299;
78 : static value_type constexpr NTSC_LUMA_GREEN = 0.587;
79 : static value_type constexpr NTSC_LUMA_BLUE = 0.114;
80 :
81 : static value_type constexpr AVERAGE_LUMA_RED = 1.0 / 3.0;
82 : static value_type constexpr AVERAGE_LUMA_GREEN = 1.0 / 3.0;
83 : static value_type constexpr AVERAGE_LUMA_BLUE = 1.0 / 3.0;
84 :
85 : class element_ref
86 : {
87 : public:
88 555378 : element_ref(matrix<T, SIZE> & m, size_type j, size_type i)
89 : : f_matrix(m)
90 : //, f_j(j)
91 : //, f_i(i)
92 555378 : , f_offset(i + j * m.columns())
93 : {
94 555378 : if(i >= m.f_columns)
95 : {
96 0 : throw std::out_of_range("used [] operator with too large a row number");
97 : }
98 555378 : }
99 :
100 221701 : element_ref & operator = (value_type const v)
101 : {
102 221701 : f_matrix.f_vector[f_offset] = v;
103 :
104 221701 : return *this;
105 : }
106 :
107 1588 : bool operator == (value_type const v) const
108 : {
109 1588 : return f_matrix.f_vector[f_offset] == v;
110 : }
111 :
112 : bool operator != (value_type const v) const
113 : {
114 : return f_matrix.f_vector[f_offset] != v;
115 : }
116 :
117 : bool operator < (value_type const v) const
118 : {
119 : return f_matrix.f_vector[f_offset] < v;
120 : }
121 :
122 : bool operator <= (value_type const v) const
123 : {
124 : return f_matrix.f_vector[f_offset] <= v;
125 : }
126 :
127 : bool operator > (value_type const v) const
128 : {
129 : return f_matrix.f_vector[f_offset] > v;
130 : }
131 :
132 : bool operator >= (value_type const v) const
133 : {
134 : return f_matrix.f_vector[f_offset] >= v;
135 : }
136 :
137 332089 : operator T () const
138 : {
139 332089 : return f_matrix.f_vector[f_offset];
140 : }
141 :
142 : private:
143 : matrix<T, SIZE> & f_matrix;
144 : size_type f_offset;
145 : };
146 :
147 : class row_ref
148 : {
149 : public:
150 555378 : row_ref(matrix<T, SIZE> & m, size_type j)
151 : : f_matrix(m)
152 555378 : , f_j(j)
153 : {
154 555378 : if(j >= m.f_rows)
155 : {
156 0 : throw std::out_of_range("used [] operator with too large a row number");
157 : }
158 555378 : }
159 :
160 555378 : element_ref operator [] (size_type i)
161 : {
162 555378 : element_ref r(f_matrix, f_j, i);
163 555378 : return r;
164 : }
165 :
166 : value_type operator [] (size_type i) const
167 : {
168 : return f_matrix.f_vector[i + f_j * f_matrix.columns()];
169 : }
170 :
171 : private:
172 : matrix<T, SIZE> & f_matrix;
173 : size_type f_j; // row
174 : };
175 :
176 : class const_row_ref
177 : {
178 : public:
179 : const_row_ref(matrix<T, SIZE> const & m, size_type j)
180 : : f_matrix(m)
181 : , f_j(j)
182 : {
183 : if(j >= m.f_rows)
184 : {
185 : throw std::out_of_range("used [] operator with too large a row number");
186 : }
187 : }
188 :
189 : value_type operator [] (size_type i) const
190 : {
191 : if(i >= f_matrix.f_columns)
192 : {
193 : throw std::out_of_range("used [] operator with too large a row number");
194 : }
195 : return f_matrix.f_vector[i + f_j * f_matrix.columns()];
196 : }
197 :
198 : private:
199 : matrix<T, SIZE> const & f_matrix;
200 : size_type f_j; // row
201 : };
202 :
203 4 : matrix<T, SIZE>()
204 4 : {
205 4 : }
206 :
207 95299 : matrix<T, SIZE>(size_type rows, size_type columns)
208 : : f_rows(rows)
209 : , f_columns(columns)
210 95299 : , f_vector(rows * columns)
211 : {
212 95299 : initialize();
213 95299 : }
214 :
215 : template<typename V, typename SZ>
216 : matrix<T, SIZE>(matrix<V, SZ> const & rhs)
217 : : f_rows(rhs.f_rows)
218 : , f_columns(rhs.f_columns)
219 : , f_vector(rhs.f_vector.size())
220 : {
221 : initialize();
222 : }
223 :
224 : template<typename V, typename SZ>
225 : matrix<T, SIZE> & operator = (matrix<V, SZ> const & rhs) //= default;
226 : {
227 : f_rows = rhs.f_rows;
228 : f_columns = rhs.f_columns;
229 : f_vector = rhs.f_vector;
230 : f_luma_red = rhs.f_luma_red;
231 : f_luma_green = rhs.f_luma_green;
232 : f_luma_blue = rhs.f_luma_blue;
233 : #ifdef _DEBUG
234 : f_last_hue_matrix = rhs.f_last_hue_matrix;
235 : #endif
236 : }
237 :
238 : template<typename V, typename SZ>
239 : matrix<T, SIZE> & operator = (matrix<V, SZ> const && rhs) //= default;
240 : {
241 : f_rows = std::move(rhs.f_rows);
242 : f_columns = std::move(rhs.f_columns);
243 : f_vector = std::move(rhs.f_vector);
244 : f_luma_red = std::move(rhs.f_luma_red);
245 : f_luma_green = std::move(rhs.f_luma_green);
246 : f_luma_blue = std::move(rhs.f_luma_blue);
247 : #ifdef _DEBUG
248 : f_last_hue_matrix = std::move(rhs.f_last_hue_matrix);
249 : #endif
250 : }
251 :
252 : bool empty() const
253 : {
254 : return f_rows == 0 || f_columns == 0;
255 : }
256 :
257 : bool is_diagonal() const
258 : {
259 : for(size_type j(0); j < f_rows; ++j)
260 : {
261 : size_type const joffset(j * f_columns);
262 : for(size_type i(0); i < f_columns; ++i)
263 : {
264 : if(i != j
265 : && f_vector[i + joffset] != static_cast<value_type>(0)) // TODO: use a compare function with epsilon instead
266 : {
267 : return false;
268 : }
269 : }
270 : }
271 : return true;
272 : }
273 :
274 72 : size_type rows() const
275 : {
276 72 : return f_rows;
277 : }
278 :
279 555450 : size_type columns() const
280 : {
281 555450 : return f_columns;
282 : }
283 :
284 3 : void swap(matrix<T, SIZE> & rhs)
285 : {
286 3 : std::swap(f_rows, rhs.f_rows);
287 3 : std::swap(f_columns, rhs.f_columns);
288 3 : f_vector.swap(rhs.f_vector);
289 :
290 : // color related parameters
291 : //
292 3 : std::swap(f_luma_red, rhs.f_luma_red);
293 3 : std::swap(f_luma_green, rhs.f_luma_green);
294 3 : std::swap(f_luma_blue, rhs.f_luma_blue);
295 : #ifdef _DEBUG
296 3 : f_last_hue_matrix.swap(rhs.f_last_hue_matrix);
297 : #endif
298 3 : }
299 :
300 95299 : void initialize()
301 : {
302 95299 : if(f_rows == f_columns)
303 : {
304 81297 : identity();
305 : }
306 : else
307 : {
308 14002 : clear();
309 : }
310 95299 : }
311 :
312 14005 : void clear()
313 : {
314 14005 : std::fill(f_vector.begin(), f_vector.end(), value_type());
315 14005 : }
316 :
317 81297 : void identity()
318 : {
319 406061 : for(size_type j(0); j < f_rows; ++j)
320 : {
321 324764 : size_type const joffset(j * f_columns);
322 1622924 : for(size_type i(0); i < f_columns; ++i)
323 : {
324 1298160 : f_vector[i + joffset] =
325 : i == j
326 1298160 : ? static_cast<value_type>(1)
327 : : value_type();
328 : }
329 : }
330 81297 : }
331 :
332 555378 : row_ref operator [] (size_type row)
333 : {
334 555378 : row_ref r(*this, row);
335 555378 : return r;
336 : }
337 :
338 : const_row_ref operator [] (size_type row) const
339 : {
340 : const_row_ref r(*this, row);
341 : return r;
342 : }
343 :
344 : template<class S>
345 3 : matrix<T, SIZE> operator * (S const & scalar) const
346 : {
347 6 : matrix<T, SIZE> t(*this);
348 6 : return t *= scalar;
349 : }
350 :
351 : // At this point I don't know how to make that work...
352 : // template<class S, typename V, typename SZ> friend
353 : // matrix<V, SZ> operator * (S const & scalar, matrix<V, SZ> const & m)
354 : // {
355 : // return m * scalar;
356 : // }
357 :
358 : template<class S>
359 4 : matrix<T, SIZE> & operator *= (S const & scalar)
360 : {
361 20 : for(size_type j(0); j < f_rows; ++j)
362 : {
363 16 : size_type const joffset(j * f_columns);
364 80 : for(size_type i(0); i < f_columns; ++i)
365 : {
366 64 : f_vector[i + joffset] *= scalar;
367 : }
368 : }
369 :
370 4 : return *this;
371 : }
372 :
373 : template<typename V, typename SZ>
374 48007 : matrix<T, SIZE> operator * (matrix<V, SZ> const & m) const
375 : {
376 48007 : if(f_columns != m.f_rows)
377 : {
378 : // enhance that support with time
379 : // (i.e. the number of columns of one matrix has to match the
380 : // number of rows of the other, they don't need to be square
381 : // matrices of the same size)
382 : //
383 0 : std::cerr << "this " << f_rows << "x" << f_columns
384 0 : << " and rhs " << m.f_rows << "x" << m.f_columns << "\n";
385 0 : throw std::runtime_error("matrices of incompatible sizes for a multiplication");
386 : }
387 :
388 : // temporary buffer
389 : //
390 48007 : matrix<T, SIZE> t(f_rows, m.f_columns);
391 :
392 240035 : for(size_type j(0); j < f_rows; ++j)
393 : {
394 192028 : size_type const joffset(j * f_columns);
395 876140 : for(size_type i(0); i < m.f_columns; ++i)
396 : {
397 684112 : value_type sum = value_type();
398 :
399 : // k goes from 0 to (f_columns == m.f_rows)
400 : //
401 3420560 : for(size_type k(0); k < m.f_rows; ++k) // sometimes it's X and sometimes it's Y
402 : {
403 5472896 : sum += f_vector[k + joffset]
404 2736448 : * m.f_vector[i + k * m.f_columns];
405 : }
406 684112 : t.f_vector[i + j * t.f_columns] = sum;
407 : }
408 : }
409 :
410 48007 : return t;
411 : }
412 :
413 : template<class V>
414 7002 : matrix<T, SIZE> & operator *= (matrix<V> const & m)
415 : {
416 7002 : return *this = *this * m;
417 : }
418 :
419 : template<class S>
420 1 : matrix<T, SIZE> operator / (S const & scalar) const
421 : {
422 2 : matrix<T, SIZE> t(*this);
423 2 : return t /= scalar;
424 : }
425 :
426 : template<class S>
427 2 : matrix<T, SIZE> & operator /= (S const & scalar)
428 : {
429 10 : for(size_type y(0); y < f_rows; ++y)
430 : {
431 40 : for(size_type x(0); x < f_columns; ++x)
432 : {
433 32 : size_type const idx(x + y * f_columns);
434 32 : f_vector[idx] /= scalar;
435 : }
436 : }
437 :
438 2 : return *this;
439 : }
440 :
441 : template<typename V, typename SZ>
442 7001 : matrix<T, SIZE> operator / (matrix<V, SZ> const & m) const
443 : {
444 : // temporary buffer
445 : //
446 14002 : matrix<T, SIZE> t(m);
447 7001 : t.inverse();
448 14002 : return *this * t;
449 : }
450 :
451 : template<typename V, typename SZ>
452 1 : matrix<T, SIZE> & operator /= (matrix<V, SZ> const & m)
453 : {
454 2 : matrix<T, SIZE> t(m);
455 1 : t.inverse();
456 2 : return *this *= t;
457 : }
458 :
459 : /** \brief Compute the inverse of `this` matrix if possible.
460 : *
461 : * This function computes the matrix determinant to see whether
462 : * it can be inverted. If so, it computes the inverse and becomes
463 : * that inverse.
464 : *
465 : * The function returns false if the inverse cannot be calculated
466 : * and the matrix remains unchanged.
467 : *
468 : * $$A^{-1} = {1 \over det(A)} adj(A)$$
469 : *
470 : * \return true if the inverse computation was successful.
471 : */
472 7003 : bool inverse()
473 : {
474 7003 : if(f_rows != 4
475 7003 : || f_columns != 4)
476 : {
477 0 : double const det(determinant());
478 : #pragma GCC diagnostic push
479 : #pragma GCC diagnostic ignored "-Wfloat-equal"
480 0 : if(det == static_cast<value_type>(0))
481 : {
482 0 : return false;
483 : }
484 : #pragma GCC diagnostic pop
485 0 : snapdev::matrix<double> adj(adjugate());
486 :
487 0 : *this = adj * (static_cast<value_type>(1) / det);
488 0 : return true;
489 : }
490 :
491 : // the following is very specific to a 4x4 matrix...
492 : //
493 7003 : value_type temp[4][8], *r[4], m[5];
494 :
495 7003 : r[0] = temp[0];
496 7003 : r[1] = temp[1];
497 7003 : r[2] = temp[2];
498 7003 : r[3] = temp[3];
499 :
500 7003 : temp[0][0] = f_vector[0 + 0 * 4];
501 7003 : temp[0][1] = f_vector[0 + 1 * 4];
502 7003 : temp[0][2] = f_vector[0 + 2 * 4];
503 7003 : temp[0][3] = f_vector[0 + 3 * 4];
504 :
505 7003 : temp[1][0] = f_vector[1 + 0 * 4];
506 7003 : temp[1][1] = f_vector[1 + 1 * 4];
507 7003 : temp[1][2] = f_vector[1 + 2 * 4];
508 7003 : temp[1][3] = f_vector[1 + 3 * 4];
509 :
510 7003 : temp[2][0] = f_vector[2 + 0 * 4];
511 7003 : temp[2][1] = f_vector[2 + 1 * 4];
512 7003 : temp[2][2] = f_vector[2 + 2 * 4];
513 7003 : temp[2][3] = f_vector[2 + 3 * 4];
514 :
515 7003 : temp[3][0] = f_vector[3 + 0 * 4];
516 7003 : temp[3][1] = f_vector[3 + 1 * 4];
517 7003 : temp[3][2] = f_vector[3 + 2 * 4];
518 7003 : temp[3][3] = f_vector[3 + 3 * 4];
519 :
520 : // can we do it?!
521 7003 : value_type abs_a = std::abs(r[3][0]);
522 7003 : value_type abs_b = std::abs(r[2][0]);
523 7003 : if(abs_a > abs_b) {
524 0 : value_type *swap = r[3];
525 0 : r[3] = r[2];
526 0 : r[2] = swap;
527 0 : abs_b = abs_a;
528 : }
529 :
530 7003 : abs_a = std::abs(r[1][0]);
531 7003 : if(abs_b > abs_a)
532 : {
533 7000 : value_type *swap = r[2];
534 7000 : r[2] = r[1];
535 7000 : r[1] = swap;
536 7000 : abs_a = abs_b;
537 : }
538 :
539 7003 : abs_b = std::abs(r[0][0]);
540 7003 : if(abs_a > abs_b)
541 : {
542 1000 : value_type *swap = r[1];
543 1000 : r[1] = r[0];
544 1000 : r[0] = swap;
545 1000 : abs_b = abs_a;
546 : }
547 :
548 : #pragma GCC diagnostic push
549 : #pragma GCC diagnostic ignored "-Wfloat-equal"
550 7003 : if(abs_b == 0.0f)
551 : {
552 0 : return false;
553 : }
554 :
555 7003 : temp[0][4] = temp[1][5] = temp[2][6] = temp[3][7] = 1.0f;
556 84036 : temp[0][5] =
557 7003 : temp[0][6] =
558 7003 : temp[0][7] =
559 7003 : temp[1][4] =
560 7003 : temp[1][6] =
561 7003 : temp[1][7] =
562 7003 : temp[2][4] =
563 7003 : temp[2][5] =
564 7003 : temp[2][7] =
565 7003 : temp[3][4] =
566 7003 : temp[3][5] =
567 7003 : temp[3][6] = 0.0f;
568 :
569 : // first elimination
570 7003 : m[4] = 1.0f / r[0][0];
571 7003 : m[1] = r[1][0] * m[4];
572 7003 : m[2] = r[2][0] * m[4];
573 7003 : m[3] = r[3][0] * m[4];
574 :
575 7003 : m[4] = r[0][1];
576 7003 : if(m[4] != 0.0)
577 : {
578 7003 : r[1][1] -= m[1] * m[4];
579 7003 : r[2][1] -= m[2] * m[4];
580 7003 : r[3][1] -= m[3] * m[4];
581 : }
582 :
583 7003 : m[4] = r[0][2];
584 7003 : if(m[4] != 0.0)
585 : {
586 7003 : r[1][2] -= m[1] * m[4];
587 7003 : r[2][2] -= m[2] * m[4];
588 7003 : r[3][2] -= m[3] * m[4];
589 : }
590 :
591 7003 : m[4] = r[0][3];
592 7003 : if(m[4] != 0.0)
593 : {
594 3 : r[1][3] -= m[1] * m[4];
595 3 : r[2][3] -= m[2] * m[4];
596 3 : r[3][3] -= m[3] * m[4];
597 : }
598 :
599 7003 : m[4] = r[0][4];
600 7003 : if(m[4] != 0.0)
601 : {
602 6003 : r[1][4] -= m[1] * m[4];
603 6003 : r[2][4] -= m[2] * m[4];
604 6003 : r[3][4] -= m[3] * m[4];
605 : }
606 :
607 7003 : m[4] = r[0][5];
608 7003 : if(m[4] != 0.0)
609 : {
610 0 : r[1][5] -= m[1] * m[4];
611 0 : r[2][5] -= m[2] * m[4];
612 0 : r[3][5] -= m[3] * m[4];
613 : }
614 :
615 7003 : m[4] = r[0][6];
616 7003 : if(m[4] != 0.0)
617 : {
618 1000 : r[1][6] -= m[1] * m[4];
619 1000 : r[2][6] -= m[2] * m[4];
620 1000 : r[3][6] -= m[3] * m[4];
621 : }
622 :
623 7003 : m[4] = r[0][7];
624 7003 : if(m[4] != 0.0)
625 : {
626 0 : r[1][7] -= m[1] * m[4];
627 0 : r[2][7] -= m[2] * m[4];
628 0 : r[3][7] -= m[3] * m[4];
629 : }
630 :
631 : // can we do it?!
632 7003 : abs_a = std::abs(r[3][1]);
633 7003 : abs_b = std::abs(r[2][1]);
634 7003 : if(abs_a > abs_b)
635 : {
636 3 : value_type *swap = r[3];
637 3 : r[3] = r[2];
638 3 : r[2] = swap;
639 3 : abs_b = abs_a;
640 : }
641 7003 : abs_a = std::abs(r[1][1]);
642 7003 : if(abs_b > abs_a)
643 : {
644 7003 : value_type *swap = r[2];
645 7003 : r[2] = r[1];
646 7003 : r[1] = swap;
647 7003 : abs_a = abs_b;
648 : }
649 :
650 7003 : if(abs_a == 0.0f)
651 : {
652 0 : return false;
653 : }
654 :
655 : // second elimination
656 7003 : m[4] = 1.0f / r[1][1];
657 :
658 7003 : m[2] = r[2][1] * m[4];
659 7003 : m[3] = r[3][1] * m[4];
660 :
661 7003 : if(m[2] != 0.0f)
662 : {
663 7003 : r[2][2] -= m[2] * r[1][2];
664 7003 : r[2][3] -= m[2] * r[1][3];
665 : }
666 :
667 7003 : if(m[3] != 0.0f)
668 : {
669 3 : r[3][2] -= m[3] * r[1][2];
670 3 : r[3][3] -= m[3] * r[1][3];
671 : }
672 :
673 7003 : m[4] = r[1][4];
674 7003 : if(m[4] != 0.0f)
675 : {
676 3 : r[2][4] -= m[2] * m[4];
677 3 : r[3][4] -= m[3] * m[4];
678 : }
679 :
680 7003 : m[4] = r[1][5];
681 7003 : if(m[4] != 0.0f)
682 : {
683 7000 : r[2][5] -= m[2] * m[4];
684 7000 : r[3][5] -= m[3] * m[4];
685 : }
686 :
687 7003 : m[4] = r[1][6];
688 7003 : if(m[4] != 0.0f)
689 : {
690 0 : r[2][6] -= m[2] * m[4];
691 0 : r[3][6] -= m[3] * m[4];
692 : }
693 :
694 7003 : m[4] = r[1][7];
695 7003 : if(m[4] != 0.0f)
696 : {
697 3 : r[2][7] -= m[2] * m[4];
698 3 : r[3][7] -= m[3] * m[4];
699 : }
700 :
701 : // can we do it?!
702 7003 : abs_a = std::abs(r[3][2]);
703 7003 : abs_b = std::abs(r[2][2]);
704 7003 : if(abs_a > abs_b)
705 : {
706 3 : value_type *swap = r[3];
707 3 : r[3] = r[2];
708 3 : r[2] = swap;
709 3 : abs_b = abs_a;
710 : }
711 :
712 7003 : if(abs_b == 0.0f)
713 : {
714 0 : return false;
715 : }
716 :
717 : // third elimination
718 7003 : m[3] = r[3][2] / r[2][2];
719 :
720 7003 : r[3][3] -= m[3] * r[2][3];
721 7003 : r[3][4] -= m[3] * r[2][4];
722 7003 : r[3][5] -= m[3] * r[2][5];
723 7003 : r[3][6] -= m[3] * r[2][6];
724 7003 : r[3][7] -= m[3] * r[2][7];
725 :
726 : // can we do it?!
727 7003 : if(r[3][3] == 0.0f) {
728 0 : return false;
729 : }
730 : #pragma GCC diagnostic pop
731 :
732 : // back substitute
733 : /* 3 */
734 7003 : m[4] = 1.0f / r[3][3];
735 7003 : r[3][4] *= m[4];
736 7003 : r[3][5] *= m[4];
737 7003 : r[3][6] *= m[4];
738 7003 : r[3][7] *= m[4];
739 :
740 : /* 2 */
741 7003 : m[2] = r[2][3];
742 7003 : m[4] = 1.0f / r[2][2];
743 7003 : r[2][4] = m[4] * (r[2][4] - r[3][4] * m[2]);
744 7003 : r[2][5] = m[4] * (r[2][5] - r[3][5] * m[2]);
745 7003 : r[2][6] = m[4] * (r[2][6] - r[3][6] * m[2]);
746 7003 : r[2][7] = m[4] * (r[2][7] - r[3][7] * m[2]);
747 :
748 7003 : m[1] = r[1][3];
749 7003 : r[1][4] -= r[3][4] * m[1];
750 7003 : r[1][5] -= r[3][5] * m[1];
751 7003 : r[1][6] -= r[3][6] * m[1];
752 7003 : r[1][7] -= r[3][7] * m[1];
753 :
754 7003 : m[0] = r[0][3];
755 7003 : r[0][4] -= r[3][4] * m[0];
756 7003 : r[0][5] -= r[3][5] * m[0];
757 7003 : r[0][6] -= r[3][6] * m[0];
758 7003 : r[0][7] -= r[3][7] * m[0];
759 :
760 : /* 1 */
761 7003 : m[1] = r[1][2];
762 7003 : m[4] = 1.0f / r[1][1];
763 7003 : r[1][4] = m[4] * (r[1][4] - r[2][4] * m[1]);
764 7003 : r[1][5] = m[4] * (r[1][5] - r[2][5] * m[1]);
765 7003 : r[1][6] = m[4] * (r[1][6] - r[2][6] * m[1]);
766 7003 : r[1][7] = m[4] * (r[1][7] - r[2][7] * m[1]);
767 :
768 7003 : m[0] = r[0][2];
769 7003 : r[0][4] -= r[2][4] * m[0];
770 7003 : r[0][5] -= r[2][5] * m[0];
771 7003 : r[0][6] -= r[2][6] * m[0];
772 7003 : r[0][7] -= r[2][7] * m[0];
773 :
774 : /* 0 */
775 7003 : m[0] = r[0][1];
776 7003 : m[4] = 1.0f / r[0][0];
777 7003 : r[0][4] = m[4] * (r[0][4] - r[1][4] * m[0]);
778 7003 : r[0][5] = m[4] * (r[0][5] - r[1][5] * m[0]);
779 7003 : r[0][6] = m[4] * (r[0][6] - r[1][6] * m[0]);
780 7003 : r[0][7] = m[4] * (r[0][7] - r[1][7] * m[0]);
781 :
782 : // save in destination (we need to put that in the back substitution)
783 7003 : f_vector[0 + 0 * 4] = r[0][4];
784 7003 : f_vector[0 + 1 * 4] = r[0][5];
785 7003 : f_vector[0 + 2 * 4] = r[0][6];
786 7003 : f_vector[0 + 3 * 4] = r[0][7];
787 :
788 7003 : f_vector[1 + 0 * 4] = r[1][4];
789 7003 : f_vector[1 + 1 * 4] = r[1][5];
790 7003 : f_vector[1 + 2 * 4] = r[1][6];
791 7003 : f_vector[1 + 3 * 4] = r[1][7];
792 :
793 7003 : f_vector[2 + 0 * 4] = r[2][4];
794 7003 : f_vector[2 + 1 * 4] = r[2][5];
795 7003 : f_vector[2 + 2 * 4] = r[2][6];
796 7003 : f_vector[2 + 3 * 4] = r[2][7];
797 :
798 7003 : f_vector[3 + 0 * 4] = r[3][4];
799 7003 : f_vector[3 + 1 * 4] = r[3][5];
800 7003 : f_vector[3 + 2 * 4] = r[3][6];
801 7003 : f_vector[3 + 3 * 4] = r[3][7];
802 :
803 7003 : return true;
804 : }
805 :
806 : /** \brief Reduce a matrix by removing one row and one column.
807 : *
808 : * This function creates a minor duplicate of this matrix with column i
809 : * and row j removed.
810 : *
811 : * The minor is denoted:
812 : *
813 : * $$M_{ij}$$
814 : *
815 : * It is a matrix built from $A$ without column `i` and row `j`.
816 : *
817 : * \note
818 : * The matrix must at least be a 2x2 matrix.
819 : *
820 : * \note
821 : * There is a "minor" macro so I named this function minor_matrix().
822 : *
823 : * \param[in] row The row to remove.
824 : * \param[in] column The column to remove.
825 : *
826 : * \return The requested minor matrix.
827 : */
828 220 : matrix<T, SIZE> minor_matrix(size_type row, size_type column) const
829 : {
830 220 : if(f_rows < 2
831 220 : || f_columns < 2)
832 : {
833 0 : throw std::runtime_error("minor_matrix() must be called with a matrix which is at least 2x2, although it does not need to be a square matrix");
834 : }
835 :
836 220 : matrix<T, SIZE> p(f_rows - 1, f_columns - 1);
837 :
838 : // we loop using p sizes
839 : // the code below ensures the correct input is retrieved
840 : //
841 : // di -- destination column
842 : // dj -- destination row
843 : // si -- source column
844 : // sj -- source row
845 : //
846 701 : for(size_type dj(0); dj < p.f_rows; ++dj)
847 : {
848 1574 : for(size_type di(0); di < p.f_columns; ++di)
849 : {
850 : // here we have 4 cases:
851 : //
852 : // a11 a12 | a13 | a14 a15
853 : // a21 a22 | a23 | a24 a25
854 : // --------+-----+--------
855 : // a31 a32 | a33 | a34 a35
856 : // --------+-----+--------
857 : // a41 a42 | a43 | a44 a45
858 : // a51 a52 | a53 | a54 a55
859 : //
860 : // assuming 'row' and 'column' are set to 3 and 3, then
861 : // the graph shows the 4 cases as the 4 corners, the
862 : // center lines will be removed so they are ignored in
863 : // the source
864 : //
865 1093 : size_type const si(di < column ? di : di + 1);
866 1093 : size_type const sj(dj < row ? dj : dj + 1);
867 :
868 1093 : p.f_vector[di + dj * p.f_columns] = f_vector[si + sj * f_columns];
869 : }
870 : }
871 :
872 220 : return p;
873 : }
874 :
875 : /** \brief Calculate the determinant of this matrix.
876 : *
877 : * This function calculates the determinant of this matrix:
878 : *
879 : * $$det(A) = \sum_{\sigma \in S_n} \Big( sgn(\sigma) \prod_{i=1}^{n} a_{i,\sigma_i} \Big)$$
880 : *
881 : * The function is implemented using a recursive set of calls. It
882 : * knows how to calculate a 2x2 matrix. All the others use recursive
883 : * calls to calculate the final result.
884 : *
885 : * Let's say you have a 3x3 matrix like this:
886 : *
887 : * \code
888 : * | a11 a12 a13 |
889 : * | a21 a22 a23 |
890 : * | a31 a32 a33 |
891 : * \endcode
892 : *
893 : * If first calculates the determinant of the 2x2 matrix:
894 : *
895 : * \code
896 : * | a22 a23 | = a22 x a33 - a23 x a32
897 : * | a32 a33 |
898 : * \endcode
899 : *
900 : * Then multiply that determinant by `a11`.
901 : *
902 : * Next it calculates the determinant of the 2x2 matrix:
903 : *
904 : * \code
905 : * | a21 a23 | = a21 x a33 - a23 x a31
906 : * | a31 a33 |
907 : * \endcode
908 : *
909 : * Then multiply that determinant by `a12` and multiply by -1.
910 : *
911 : * Finally, it calculates the last 2x2 matrix at the bottom left corner.
912 : *
913 : * \code
914 : * | a21 a22 | = a21 x a32 - a22 x a31
915 : * | a31 a32 |
916 : * \endcode
917 : *
918 : * Then multiply that last determinant by `a13`.
919 : *
920 : * Finally we sum all three results and that's our determinant for a
921 : * 3x3 matrix.
922 : *
923 : * The determinant of a 4x4 matrix will be calculated in a similar
924 : * way by also calculating the determinant of all the 3x3 matrices
925 : * defined under the first row.
926 : *
927 : * Source: https://en.wikipedia.org/wiki/Determinant
928 : *
929 : * \exception runtime_error
930 : * If the matrix is not a square matrix, raise a runtime_error exception.
931 : *
932 : * \return The determinant value.
933 : */
934 229 : value_type determinant() const
935 : {
936 229 : if(f_rows != f_columns)
937 : {
938 0 : throw std::runtime_error("determinant can only be calculated for square matrices");
939 : }
940 :
941 229 : if(f_columns == 1)
942 : {
943 4 : return f_vector[0];
944 : }
945 :
946 225 : if(f_columns == 2)
947 : {
948 172 : return f_vector[0 + 0 * 2] * f_vector[1 + 1 * 2]
949 172 : - f_vector[1 + 0 * 2] * f_vector[0 + 1 * 2];
950 : }
951 :
952 53 : value_type determinant = value_type();
953 :
954 53 : value_type sign = static_cast<value_type>(1);
955 214 : for(size_type c(0); c < f_columns; ++c)
956 : {
957 : // create a minor matrix
958 : //
959 322 : matrix<T, SIZE> p(minor_matrix(0, c));
960 :
961 : // add to the determinant for that column
962 : // (number of row 0 column 'c')
963 : //
964 161 : determinant += sign
965 161 : * f_vector[c + 0 * f_columns]
966 161 : * p.determinant();
967 :
968 : // swap the sign
969 : //
970 161 : sign *= static_cast<value_type>(-1);
971 : }
972 :
973 53 : return determinant;
974 : }
975 :
976 : /** \brief Swap the rows and columns of a matrix.
977 : *
978 : * This function returns the transpose of this matrix.
979 : *
980 : * Generally noted as:
981 : *
982 : * $$A^T$$
983 : *
984 : * The definition of the transpose is:
985 : *
986 : * $$[A^T]_{ij} = [A]_{ji}$$
987 : *
988 : * The resulting matrix has a number of columns equal to 'this' matrix
989 : * rows and vice versa.
990 : *
991 : * \return A new matrix representing the transpose of 'this' matrix.
992 : */
993 6 : matrix<T, SIZE> transpose() const
994 : {
995 : // 'm' has its number of rows and columns swapped compared
996 : // to 'this'
997 6 : matrix<T, SIZE> m(f_columns, f_rows);
998 :
999 29 : for(size_type j(0); j < f_rows; ++j)
1000 : {
1001 96 : for(size_type i(0); i < f_columns; ++i)
1002 : {
1003 : // we could also have used "j + i * f_rows" on the left
1004 : // but I think it's more confusing
1005 : //
1006 73 : m.f_vector[j + i * m.f_columns] = f_vector[i + j * f_columns];
1007 : }
1008 : }
1009 :
1010 6 : return m;
1011 : }
1012 :
1013 : /** \brief This function calculates the adjugate of this matrix.
1014 : *
1015 : * \return The adjugate of this matrix.
1016 : */
1017 4 : matrix<T, SIZE> adjugate() const
1018 : {
1019 4 : if(f_rows != f_columns)
1020 : {
1021 : // is that true?
1022 : //
1023 0 : throw std::runtime_error("adjugate can only be calculated for square matrices");
1024 : }
1025 :
1026 8 : matrix<T, SIZE> r(f_rows, f_columns);
1027 :
1028 : // det(A) when A is 1x1 equals | 1 |
1029 : // which is the default 'r'
1030 : //
1031 4 : if(f_columns != 1)
1032 : {
1033 : //if(f_columns == 2)
1034 : //{
1035 : // // the 2x2 matrix is handled as a special case just so it goes
1036 : // // faster but not so much more than that
1037 : // //
1038 : // // adj | a b | = | d -b |
1039 : // // | c d | | -c a |
1040 : // //
1041 : // r.f_vector[0 + 0 * 2] = f_vector[1 + 1 * 2];
1042 : // r.f_vector[1 + 0 * 2] = -f_vector[1 + 0 * 2];
1043 : // r.f_vector[0 + 1 * 2] = -f_vector[0 + 1 * 2];
1044 : // r.f_vector[1 + 1 * 2] = f_vector[0 + 0 * 2];
1045 : //}
1046 : //else
1047 : {
1048 : // for larger matrices we use a loop and calculate the determinant
1049 : // for each new value with the "rest" of the matrix at that point
1050 : //
1051 17 : for(size_type j(0); j < f_rows; ++j)
1052 : {
1053 58 : for(size_type i(0); i < f_columns; ++i)
1054 : {
1055 90 : matrix<T, SIZE> p(minor_matrix(j, i));
1056 45 : r.f_vector[i + j * r.f_columns] = static_cast<value_type>(((i + j) & 1) == 0 ? 1 : -1) * p.determinant();
1057 : }
1058 : }
1059 :
1060 4 : return r.transpose();
1061 : }
1062 : }
1063 :
1064 0 : return r;
1065 : }
1066 :
1067 : /** \brief Add a scale to 'this' matrix.
1068 : *
1069 : * This function adds the specified scalar to the matrix. This adds
1070 : * the specified amount to all the elements of the matrix.
1071 : *
1072 : * $$[A]_{ij} \leftarrow [A]_{ij} + scalar$$
1073 : *
1074 : * \param[in] scalar The scalar to add to this matrix.
1075 : */
1076 : template<class S>
1077 1 : matrix<T, SIZE> operator + (S const & scalar) const
1078 : {
1079 2 : matrix<T, SIZE> t(*this);
1080 2 : return t += scalar;
1081 : }
1082 :
1083 : template<class S>
1084 2 : matrix<T, SIZE> & operator += (S const & scalar)
1085 : {
1086 10 : for(size_type y(0); y < f_rows; ++y)
1087 : {
1088 40 : for(size_type x(0); x < f_columns; ++x)
1089 : {
1090 32 : size_type const idx(x + y * f_columns);
1091 32 : f_vector[idx] += scalar;
1092 : }
1093 : }
1094 :
1095 2 : return *this;
1096 : }
1097 :
1098 : template<typename V, typename SZ>
1099 1 : matrix<T, SIZE> operator + (matrix<V, SZ> const & m) const
1100 : {
1101 2 : matrix<T, SIZE> t(*this);
1102 2 : return t += m;
1103 : }
1104 :
1105 : template<typename V, typename SZ>
1106 2 : matrix<T, SIZE> & operator += (matrix<V, SZ> const & m)
1107 : {
1108 2 : if(f_rows != m.f_rows
1109 2 : || f_columns != m.f_columns)
1110 : {
1111 0 : throw std::runtime_error("matrices of incompatible sizes for an addition");
1112 : }
1113 :
1114 10 : for(size_type y(0); y < f_rows; ++y)
1115 : {
1116 40 : for(size_type x(0); x < f_columns; ++x)
1117 : {
1118 32 : size_type const idx(x + y * f_columns);
1119 32 : f_vector[idx] += m.f_vector[idx];
1120 : }
1121 : }
1122 :
1123 2 : return *this;
1124 : }
1125 :
1126 : template<class S>
1127 1 : matrix<T, SIZE> operator - (S const & scalar) const
1128 : {
1129 2 : matrix<T, SIZE> t(*this);
1130 2 : return t -= scalar;
1131 : }
1132 :
1133 : template<class S>
1134 2 : matrix<T, SIZE> & operator -= (S const & scalar)
1135 : {
1136 4 : matrix<T, SIZE> t(f_rows, f_columns);
1137 :
1138 10 : for(size_type y(0); y < f_rows; ++y)
1139 : {
1140 40 : for(size_type x(0); x < f_columns; ++x)
1141 : {
1142 32 : size_type const idx(x + y * f_columns);
1143 32 : f_vector[idx] -= scalar;
1144 : }
1145 : }
1146 :
1147 4 : return *this;
1148 : }
1149 :
1150 : template<typename V, typename SZ>
1151 1 : matrix<T, SIZE> operator - (matrix<V, SZ> const & m) const
1152 : {
1153 2 : matrix<T, SIZE> t(*this);
1154 2 : return t -= m;
1155 : }
1156 :
1157 : template<class V, typename SZ>
1158 2 : matrix<T, SIZE> & operator -= (matrix<V, SZ> const & m)
1159 : {
1160 2 : if(f_rows != m.f_rows
1161 2 : || f_columns != m.f_columns)
1162 : {
1163 0 : throw std::runtime_error("matrices of incompatible sizes for a subtraction");
1164 : }
1165 :
1166 10 : for(size_type y(0); y < f_rows; ++y)
1167 : {
1168 40 : for(size_type x(0); x < f_columns; ++x)
1169 : {
1170 32 : size_type const idx(x + y * f_columns);
1171 32 : f_vector[idx] -= m.f_vector[idx];
1172 : }
1173 : }
1174 :
1175 2 : return *this;
1176 : }
1177 :
1178 :
1179 : static matrix<T, SIZE> rotation_matrix_4x4_x(double angle)
1180 : {
1181 : matrix<T, SIZE> result(4, 4);
1182 :
1183 : double const c(cos(angle));
1184 : double const s(sin(angle));
1185 :
1186 : result[1][1] = c;
1187 : result[1][2] = -s;
1188 : result[2][1] = s;
1189 : result[2][2] = c;
1190 :
1191 : return result;
1192 : }
1193 :
1194 :
1195 : static matrix<T, SIZE> rotation_matrix_4x4_y(double angle)
1196 : {
1197 : matrix<T, SIZE> result(4, 4);
1198 :
1199 : double const c(cos(angle));
1200 : double const s(sin(angle));
1201 :
1202 : result[0][0] = c;
1203 : result[0][2] = s;
1204 : result[2][0] = -s;
1205 : result[2][2] = c;
1206 :
1207 : return result;
1208 : }
1209 :
1210 :
1211 : static matrix<T, SIZE> rotation_matrix_4x4_z(double angle)
1212 : {
1213 : matrix<T, SIZE> result(4, 4);
1214 :
1215 : double const c(cos(angle));
1216 : double const s(sin(angle));
1217 :
1218 : result[0][0] = c;
1219 : result[0][1] = -s;
1220 : result[2][0] = s;
1221 : result[2][1] = c;
1222 :
1223 : return result;
1224 : }
1225 :
1226 :
1227 : /** \brief Apply a scaling factor to this matrix.
1228 : *
1229 : * This function multiplies 'this' matrix by a scaling factor and
1230 : * returns the result. 'this' does not get changed.
1231 : *
1232 : * The scale matrix looks like:
1233 : *
1234 : * $$
1235 : * \begin{bmatrix}
1236 : * brightness_r & 0 & 0 & 0
1237 : * \\ 0 & brightness_g & 0 & 0
1238 : * \\ 0 & 0 & brightness_b & 0
1239 : * \\ 0 & 0 & 0 & 1
1240 : * \end{bmatrix}
1241 : * $$
1242 : *
1243 : * The 'r', 'g', 'b' indices represent the RGB colors if one wants to
1244 : * scale just one color at a time although this function only offers
1245 : * to set all of the fields to the same value \p s.
1246 : *
1247 : * See:
1248 : * http://www.graficaobscura.com/matrix/index.html
1249 : *
1250 : * \param[in] b The scaling factor to apply to this matrix.
1251 : *
1252 : * \return 'this' x 'brightness'
1253 : */
1254 : template<typename S>
1255 1000 : matrix<T, SIZE> brightness(S const b) const
1256 : {
1257 1000 : if(f_rows != 4
1258 1000 : || f_columns != 4)
1259 : {
1260 0 : throw std::runtime_error("brightness() is only for 4x4 matrices at this time");
1261 : }
1262 :
1263 2000 : matrix<T, SIZE> m(4, 4);
1264 1000 : m[0][0] = b;
1265 1000 : m[1][1] = b;
1266 1000 : m[2][2] = b;
1267 :
1268 2000 : return *this * m;
1269 : }
1270 :
1271 :
1272 : /** \brief Apply an RGB color saturation to this matrix.
1273 : *
1274 : * This function applies the saturation matrix defined below to
1275 : * 'this' matrix. When the saturation parameter is set to zero (0)
1276 : * the function transforms all the colors to gray. When the saturation
1277 : * is set to one (1), the color does not get changed.
1278 : *
1279 : * A saturation parameter outside of the [0 .. 1] range will have
1280 : * unexpected results.
1281 : *
1282 : * The saturation matrix:
1283 : *
1284 : * $$
1285 : * \begin{bmatrix}
1286 : * L_r \times ( 1 - saturation) + saturation & L_r \times ( 1 - saturation) & L_r \times ( 1 - saturation) & 0
1287 : * \\ L_g \times ( 1 - saturation) & L_g \times ( 1 - saturation) + saturation & L_g \times ( 1 - saturation) & 0
1288 : * \\ L_b \times ( 1 - saturation) & L_b \times ( 1 - saturation) & L_b \times ( 1 - saturation) + saturation & 0
1289 : * \\ 0 & 0 & 0 & 1
1290 : * \end{bmatrix}
1291 : * $$
1292 : *
1293 : * The weights used here come from the luma matrix. See the
1294 : * get_luma_vector() function.
1295 : *
1296 : * See:
1297 : * http://www.graficaobscura.com/matrix/index.html
1298 : *
1299 : * \param[in] s How close or far to correct the color toward saturation.
1300 : *
1301 : * \return 'this' x 'scale'
1302 : */
1303 : template<typename S>
1304 6000 : matrix<T, SIZE> saturation(S const s) const
1305 : {
1306 6000 : if(f_rows != 4
1307 6000 : || f_columns != 4)
1308 : {
1309 0 : throw std::runtime_error("saturation() is only for 4x4 matrices at this time");
1310 : }
1311 :
1312 12000 : matrix<T, SIZE> m(4, 4);
1313 :
1314 6000 : value_type const ns(static_cast<value_type>(s));
1315 6000 : value_type const os(static_cast<value_type>(1) - static_cast<value_type>(s));
1316 :
1317 6000 : m[0][0] = f_luma_red * os + ns;
1318 6000 : m[0][1] = f_luma_red * os;
1319 6000 : m[0][2] = f_luma_red * os;
1320 :
1321 6000 : m[1][0] = f_luma_green * os;
1322 6000 : m[1][1] = f_luma_green * os + ns;
1323 6000 : m[1][2] = f_luma_green * os;
1324 :
1325 6000 : m[2][0] = f_luma_blue * os;
1326 6000 : m[2][1] = f_luma_blue * os;
1327 6000 : m[2][2] = f_luma_blue * os + ns;
1328 :
1329 12000 : return *this * m;
1330 : }
1331 :
1332 :
1333 : /** \brief Apply a hue correction.
1334 : *
1335 : * This function applies a hue correction to 'this' matrix.
1336 : *
1337 : * The hue correction makes use of the rotation around the red
1338 : * and green axis, then a skew to take luminance in account.
1339 : * At that point it rotates the color. Finally, the function
1340 : * reverses the skew and rotate back around the green and red
1341 : * axis.
1342 : *
1343 : * The following is the list of matrices used to rotate the hue:
1344 : *
1345 : * Rotation around the Red axis (Rr):
1346 : *
1347 : * $$
1348 : * R_r =
1349 : * \begin{bmatrix}
1350 : * 1 & 0 & 0 & 0
1351 : * \\ 0 & {1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
1352 : * \\ 0 & -{1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
1353 : * \\ 0 & 0 & 0 & 1
1354 : * \end{bmatrix}
1355 : * $$
1356 : *
1357 : * \note
1358 : * The $1 \over \sqrt 2$ is $sin ( \pi \over 4 )$ and
1359 : * $cos ( \pi \over 4 )$.
1360 : *
1361 : * Rotation around the Green axis (Rg):
1362 : *
1363 : * $$
1364 : * R_g =
1365 : * \begin{bmatrix}
1366 : * {\sqrt 2 \over \sqrt 3} & 0 & {1 \over \sqrt 3} & 0
1367 : * \\ 0 & 1 & 0 & 0
1368 : * \\ -{1 \over \sqrt 3} & 0 & {\sqrt 2 \over \sqrt 3} & 0
1369 : * \\ 0 & 0 & 0 & 1
1370 : * \end{bmatrix}
1371 : * $$
1372 : *
1373 : * \note
1374 : * The $\sqrt 2 \over \sqrt 3$ and $1 \over \sqrt 3$ are sin and
1375 : * cos as well. These take the first rotation in account (i.e.
1376 : * so it is again a 45° angle but applied after the 45° angle
1377 : * applied to the red axis.)
1378 : *
1379 : * Combine both rotations:
1380 : *
1381 : * $$
1382 : * R_{rg} = R_r R_g
1383 : * $$
1384 : *
1385 : * Since colors are linear but not at a similar angle, we want to
1386 : * unskew them for which we need to use the luminance. So here we
1387 : * compute the luminance of the color.
1388 : *
1389 : * $$
1390 : * \begin{bmatrix}
1391 : * L_r
1392 : * \\ L_g
1393 : * \\ L_b
1394 : * \\ 0
1395 : * \end{bmatrix}
1396 : * =
1397 : * R_{rg}
1398 : * \begin{bmatrix}
1399 : * W_r
1400 : * \\ W_g
1401 : * \\ W_b
1402 : * \\ 0
1403 : * \end{bmatrix}
1404 : * $$
1405 : *
1406 : * Now we define a rotation matrix for the blue axis. This one uses
1407 : * a variable as the angle. This is the actual rotation angle which
1408 : * can go from 0 to $2 \pi$:
1409 : *
1410 : * $$
1411 : * R_b =
1412 : * \begin{bmatrix}
1413 : * cos( \theta ) & sin( \theta ) & 0 & 0
1414 : * \\ -sin( \theta ) & cos( \theta ) & 0 & 0
1415 : * \\ 0 & 0 & 1 & 0
1416 : * \\ 0 & 0 & 0 & 1
1417 : * \end{bmatrix}
1418 : * $$
1419 : *
1420 : * Now we have all the matrices to calculate the hue rotation
1421 : * of all the components of an image:
1422 : *
1423 : * $$
1424 : * H = R_{rg} S R_b S^{-1} R_{rg}^{-1}
1425 : * $$
1426 : *
1427 : * $H$ can then be used as in:
1428 : *
1429 : * $$
1430 : * \begin{bmatrix}
1431 : * R'
1432 : * \\ G'
1433 : * \\ B'
1434 : * \end{bmatrix}
1435 : * =
1436 : * H
1437 : * \begin{bmatrix}
1438 : * R
1439 : * \\ G
1440 : * \\ B
1441 : * \end{bmatrix}
1442 : * $$
1443 : *
1444 : * See:
1445 : * http://www.graficaobscura.com/matrix/index.html
1446 : *
1447 : * We can also rewrite the hue matrix as follow:
1448 : *
1449 : * $$
1450 : * H = cos ( \theta ) C + sin ( \theta ) S + T
1451 : * $$
1452 : *
1453 : * Where C is the cosine matrix, S is the sine matrix and T is an
1454 : * additional translation.
1455 : *
1456 : * I found an example on Microsoft (see link below) which I suspect
1457 : * does not include the color skewing. In other word, it completely
1458 : * ignores the color luminance weights which gives invalid results
1459 : * (albeit in most cases relatively good.)
1460 : *
1461 : * $$
1462 : * H =
1463 : * cos ( \theta )
1464 : * \begin{bmatrix}
1465 : * 0.787 & -0.213 & -0.213
1466 : * \\ -0.715 & 0.285 & -0.715
1467 : * \\ -0.072 & -0.072 & 0.928
1468 : * \end{bmatrix}
1469 : * +
1470 : * sin ( \theta )
1471 : * \begin{bmatrix}
1472 : * -0.213 & 0.143 & -0.787
1473 : * \\ -0.715 & 0.140 & 0.715
1474 : * \\ 0.928 & -0.283 & 0.072
1475 : * \end{bmatrix}
1476 : * +
1477 : * \begin{bmatrix}
1478 : * 0.213 & 0.213 & 0.213
1479 : * \\ 0.715 & 0.715 & 0.715
1480 : * \\ 0.072 & 0.072 & 0.072
1481 : * \end{bmatrix}
1482 : * $$
1483 : *
1484 : * The correct version is the one we actually use and it goes like this
1485 : * with high precision (double floating points) and 4x4 matrices:
1486 : *
1487 : * IMPORTANT: the weights change depending on the selected luma.
1488 : * If the user makes use of a luma which is not one supported by
1489 : * default, the algorithm falls back to the fully dynamic version
1490 : * of the hue computation.
1491 : *
1492 : * $$
1493 : * H =
1494 : * cos ( \theta )
1495 : * \begin{bmatrix}
1496 : * 0.943571345820976240 & -0.056428654178995265 & -0.056428654178995265 & 0
1497 : * \\ -0.189552583569840000 & 0.810447416430107000 & -0.189552583569853000 & 0
1498 : * \\ -0.754018762251120000 & -0.754018762251113000 & 0.245981237748847000 & 0
1499 : * \\ 0 & 0 & 0 & 0
1500 : * \end{bmatrix}
1501 : * +
1502 : * sin ( \theta )
1503 : * \begin{bmatrix}
1504 : * 0.32589470021007605 & 0.90324496939967125 & -0.25145556897954369 & 0
1505 : * \\ -0.98010410586906000 & -0.40275383667945400 & 0.17459643251014600 & 0
1506 : * \\ 0.65420940565900000 & -0.50049113272020600 & 0.07685913646939400 & 0
1507 : * \\ 0 & 0 & 0 & 0
1508 : * \end{bmatrix}
1509 : * +
1510 : * \begin{bmatrix}
1511 : * 0.056428654178995 & 0.056428654178995 & 0.056428654178995 & 0
1512 : * \\ 0.189552583569860 & 0.189552583569860 & 0.189552583569860 & 0
1513 : * \\ 0.754018762251160 & 0.754018762251160 & 0.754018762251160 & 0
1514 : * \\ 0 & 0 & 0 & 1
1515 : * \end{bmatrix}
1516 : * $$
1517 : *
1518 : * This is exactly the matrix we use below, although we directly apply the
1519 : * matrix additions and multiplications to generate the matrix as quickly
1520 : * as possible (we calculate just the necessary component and avoid many
1521 : * copies so it ought to be much faster.)
1522 : *
1523 : * \note
1524 : * It looks like the Microsoft matrices are reversed, so affecting BGR
1525 : * instead of RGB. Our matrices expects RGBX, but we adjust the
1526 : * multiplication as required to handle BGR, XBGR, XRGB, etc. We will
1527 : * change our matrix if required by OpenGL but as far as I know it's
1528 : * already defined in the correct direction.
1529 : *
1530 : * See:
1531 : * https://msdn.microsoft.com/en-us/library/windows/desktop/hh706342(v=vs.85).aspx
1532 : *
1533 : * \note
1534 : * To verify that the angle is correct, one can use this Wikipedia
1535 : * page. As we can see, from Red, you add 30° to get yellow, 120°
1536 : * to get green, etc. Obviously you can go also backward with
1537 : * negative angles.
1538 : * https://en.wikipedia.org/wiki/Hue
1539 : *
1540 : * \note
1541 : * For test purposes, which version of the matrix is used is defined
1542 : * by a variable that can be retrieved using the get_last_hue_matrix().
1543 : * The function is only available in debug mode.
1544 : *
1545 : * \param[in] h The amount of rotation, an angle in radian.
1546 : *
1547 : * \return 'this' rotated around the hue axis by \p h
1548 : */
1549 : template<typename S>
1550 6000 : matrix<T, SIZE> hue(S const h) const
1551 : {
1552 6000 : if(f_rows != 4
1553 6000 : || f_columns != 4)
1554 : {
1555 0 : throw std::runtime_error("hue() is only for 4x4 matrices at this time");
1556 : }
1557 :
1558 6000 : value_type const rot_cos(std::cos(static_cast<value_type>(h)));
1559 6000 : value_type const rot_sin(std::sin(static_cast<value_type>(h)));
1560 :
1561 6000 : if(std::fabs(f_luma_red - HDTV_LUMA_RED ) < 0.0001
1562 1000 : && std::fabs(f_luma_green - HDTV_LUMA_GREEN) < 0.0001
1563 1000 : && std::fabs(f_luma_blue - HDTV_LUMA_BLUE ) < 0.0001)
1564 : {
1565 : // this matrix makes use of the HDTV luma
1566 : //
1567 2000 : matrix<T, SIZE> hue_matrix(4, 4);
1568 :
1569 1000 : hue_matrix[0][0] = 0.85089741314769186 * rot_cos + 0.39419567713872435 * rot_sin + 0.14910258685230815;
1570 1000 : hue_matrix[0][1] = -0.14910258685230816 * rot_cos + 0.97154594632835023 * rot_sin + 0.14910258685230815;
1571 1000 : hue_matrix[0][2] = -0.14910258685230816 * rot_cos - 0.18315459205090151 * rot_sin + 0.14910258685230815;
1572 :
1573 1000 : hue_matrix[1][0] = -0.08406523610970199 * rot_cos - 0.93399661436972614 * rot_sin + 0.08406523610970204;
1574 1000 : hue_matrix[1][1] = 0.91593476389029794 * rot_cos - 0.35664634518010058 * rot_sin + 0.08406523610970204;
1575 1000 : hue_matrix[1][2] = -0.08406523610970201 * rot_cos + 0.22070392400952521 * rot_sin + 0.08406523610970204;
1576 :
1577 1000 : hue_matrix[2][0] = -0.76683217703798975 * rot_cos + 0.53980093723100184 * rot_sin + 0.76683217703798978;
1578 1000 : hue_matrix[2][1] = -0.76683217703798980 * rot_cos - 0.61489960114824952 * rot_sin + 0.76683217703798978;
1579 1000 : hue_matrix[2][2] = 0.23316782296201015 * rot_cos - 0.03754933195862373 * rot_sin + 0.76683217703798978;
1580 :
1581 : #ifdef _DEBUG
1582 1000 : f_last_hue_matrix = "hdtv";
1583 : #endif
1584 :
1585 1000 : return *this * hue_matrix;
1586 : }
1587 :
1588 5000 : if(std::fabs(f_luma_red - LED_LUMA_RED ) < 0.0001
1589 1000 : && std::fabs(f_luma_green - LED_LUMA_GREEN) < 0.0001
1590 1000 : && std::fabs(f_luma_blue - LED_LUMA_BLUE ) < 0.0001)
1591 : {
1592 : // this matrix make use of the LED luma
1593 : //
1594 2000 : matrix<T, SIZE> hue_matrix(4, 4);
1595 :
1596 1000 : hue_matrix[0][0] = 0.86455583487454547 * rot_cos + 0.40703991394281032 * rot_sin + 0.13544416512545457;
1597 1000 : hue_matrix[0][1] = -0.13544416512545459 * rot_cos + 0.98439018313243625 * rot_sin + 0.13544416512545460;
1598 1000 : hue_matrix[0][2] = -0.13544416512545459 * rot_cos - 0.17031035524681553 * rot_sin + 0.13544416512545460;
1599 :
1600 1000 : hue_matrix[1][0] = -0.07977101160856729 * rot_cos - 0.95224727296282565 * rot_sin + 0.07977101160856727;
1601 1000 : hue_matrix[1][1] = 0.92022898839143270 * rot_cos - 0.37489700377320009 * rot_sin + 0.07977101160856730;
1602 1000 : hue_matrix[1][2] = -0.07977101160856729 * rot_cos + 0.20245326541642572 * rot_sin + 0.07977101160856730;
1603 :
1604 1000 : hue_matrix[2][0] = -0.78478482326597805 * rot_cos + 0.54520735902001539 * rot_sin + 0.78478482326597820;
1605 1000 : hue_matrix[2][1] = -0.78478482326597812 * rot_cos - 0.60949317935923603 * rot_sin + 0.78478482326597832;
1606 1000 : hue_matrix[2][2] = 0.21521517673402187 * rot_cos - 0.03214291016961022 * rot_sin + 0.78478482326597832;
1607 :
1608 : #ifdef _DEBUG
1609 1000 : f_last_hue_matrix = "led";
1610 : #endif
1611 :
1612 1000 : return *this * hue_matrix;
1613 : }
1614 :
1615 4000 : if(std::fabs(f_luma_red - CRT_LUMA_RED ) < 0.0001
1616 1000 : && std::fabs(f_luma_green - CRT_LUMA_GREEN) < 0.0001
1617 1000 : && std::fabs(f_luma_blue - CRT_LUMA_BLUE ) < 0.0001)
1618 : {
1619 : // this matrix make use of the CRT luma
1620 : //
1621 2000 : matrix<T, SIZE> hue_matrix(4, 4);
1622 :
1623 1000 : hue_matrix[0][0] = 0.943571345820976240 * rot_cos + 0.32589470021007605 * rot_sin + 0.056428654178995;
1624 1000 : hue_matrix[0][1] = -0.056428654178995265 * rot_cos + 0.90324496939967125 * rot_sin + 0.056428654178995;
1625 1000 : hue_matrix[0][2] = -0.056428654178995265 * rot_cos - 0.25145556897954369 * rot_sin + 0.056428654178995;
1626 :
1627 1000 : hue_matrix[1][0] = -0.189552583569840000 * rot_cos - 0.98010410586906000 * rot_sin + 0.189552583569860;
1628 1000 : hue_matrix[1][1] = 0.810447416430107000 * rot_cos - 0.40275383667945400 * rot_sin + 0.189552583569860;
1629 1000 : hue_matrix[1][2] = -0.189552583569853000 * rot_cos + 0.17459643251014600 * rot_sin + 0.189552583569860;
1630 :
1631 1000 : hue_matrix[2][0] = -0.754018762251120000 * rot_cos + 0.65420940565900000 * rot_sin + 0.754018762251160;
1632 1000 : hue_matrix[2][1] = -0.754018762251113000 * rot_cos - 0.50049113272020600 * rot_sin + 0.754018762251160;
1633 1000 : hue_matrix[2][2] = 0.245981237748847000 * rot_cos + 0.07685913646939400 * rot_sin + 0.754018762251160;
1634 :
1635 : #ifdef _DEBUG
1636 1000 : f_last_hue_matrix = "crt";
1637 : #endif
1638 :
1639 1000 : return *this * hue_matrix;
1640 : }
1641 :
1642 3000 : if(std::fabs(f_luma_red - NTSC_LUMA_RED ) < 0.0001
1643 1000 : && std::fabs(f_luma_green - NTSC_LUMA_GREEN) < 0.0001
1644 1000 : && std::fabs(f_luma_blue - NTSC_LUMA_BLUE ) < 0.0001)
1645 : {
1646 : // this matrix make use of the NTSC luma
1647 : //
1648 2000 : matrix<T, SIZE> hue_matrix(4, 4);
1649 :
1650 1000 : hue_matrix[0][0] = 0.97667266520552899 * rot_cos + 0.35888772800180165 * rot_sin + 0.02332733479447109;
1651 1000 : hue_matrix[0][1] = -0.02332733479447109 * rot_cos + 0.93623799719142759 * rot_sin + 0.02332733479447109;
1652 1000 : hue_matrix[0][2] = -0.02332733479447108 * rot_cos - 0.21846254118782418 * rot_sin + 0.02332733479447109;
1653 :
1654 1000 : hue_matrix[1][0] = -0.17753044304672443 * rot_cos - 1.02526720325074270 * rot_sin + 0.17753044304672438;
1655 1000 : hue_matrix[1][1] = 0.82246955695327556 * rot_cos - 0.44791693406111712 * rot_sin + 0.17753044304672441;
1656 1000 : hue_matrix[1][2] = -0.17753044304672441 * rot_cos + 0.12943333512850867 * rot_sin + 0.17753044304672441;
1657 :
1658 1000 : hue_matrix[2][0] = -0.79914222215880441 * rot_cos + 0.66637947524894110 * rot_sin + 0.79914222215880448;
1659 1000 : hue_matrix[2][1] = -0.79914222215880448 * rot_cos - 0.48832106313031034 * rot_sin + 0.79914222215880459;
1660 1000 : hue_matrix[2][2] = 0.20085777784119549 * rot_cos + 0.08902920605931547 * rot_sin + 0.79914222215880459;
1661 :
1662 : #ifdef _DEBUG
1663 1000 : f_last_hue_matrix = "ntsc";
1664 : #endif
1665 :
1666 1000 : return *this * hue_matrix;
1667 : }
1668 :
1669 2000 : if(std::fabs(f_luma_red - AVERAGE_LUMA_RED ) < 0.0001
1670 1000 : && std::fabs(f_luma_green - AVERAGE_LUMA_GREEN) < 0.0001
1671 1000 : && std::fabs(f_luma_blue - AVERAGE_LUMA_BLUE ) < 0.0001)
1672 : {
1673 : // this matrix make use of the average luma
1674 : // (in other words, it makes no luma correction at all)
1675 : //
1676 2000 : matrix<T, SIZE> hue_matrix(4, 4);
1677 :
1678 1000 : hue_matrix[0][0] = 1.88796748671567113 * rot_cos + 0.76774179094706859 * rot_sin - 0.88796748671567094;
1679 1000 : hue_matrix[0][1] = 0.88796748671567144 * rot_cos + 1.34509206013669466 * rot_sin - 0.88796748671567149;
1680 1000 : hue_matrix[0][2] = 0.88796748671567144 * rot_cos + 0.19039152175744295 * rot_sin - 0.88796748671567149;
1681 :
1682 1000 : hue_matrix[1][0] = -0.27909984885071244 * rot_cos - 2.01889870048836475 * rot_sin + 0.27909984885071226;
1683 1000 : hue_matrix[1][1] = 0.72090015114928749 * rot_cos - 1.44154843129873957 * rot_sin + 0.27909984885071243;
1684 1000 : hue_matrix[1][2] = -0.27909984885071243 * rot_cos - 0.86419816210911379 * rot_sin + 0.27909984885071243;
1685 :
1686 1000 : hue_matrix[2][0] = -1.60886763786495844 * rot_cos + 1.25115690954129627 * rot_sin + 1.60886763786495757;
1687 1000 : hue_matrix[2][1] = -1.60886763786495881 * rot_cos + 0.09645637116204509 * rot_sin + 1.60886763786495868;
1688 1000 : hue_matrix[2][2] = -0.60886763786495889 * rot_cos + 0.67380664035167087 * rot_sin + 1.60886763786495868;
1689 :
1690 : #ifdef _DEBUG
1691 1000 : f_last_hue_matrix = "average";
1692 : #endif
1693 :
1694 1000 : return *this * hue_matrix;
1695 : }
1696 :
1697 : #ifdef _DEBUG
1698 1000 : f_last_hue_matrix = "dynamic";
1699 : #endif
1700 :
1701 : {
1702 :
1703 : // the full computation, it works, it's just slower than using
1704 : // a pre-calculated matrix
1705 :
1706 : // $R_r$ -- rotation around red axis (inverse rotation around X axis)
1707 : //
1708 2000 : matrix<T, SIZE> r_r(4, 4);
1709 1000 : value_type const inv_sqrt_2 = static_cast<value_type>(1) / std::sqrt(static_cast<value_type>(2));
1710 1000 : r_r[1][1] = inv_sqrt_2;
1711 1000 : r_r[1][2] = inv_sqrt_2;
1712 1000 : r_r[2][1] = -inv_sqrt_2;
1713 1000 : r_r[2][2] = inv_sqrt_2;
1714 :
1715 : //std::cerr << "R_r = " << r_r << "\n";
1716 :
1717 : // $R_g$ -- rotation around green axis (inverse rotation around Y axis)
1718 : //
1719 2000 : matrix<T, SIZE> r_g(4, 4);
1720 1000 : value_type const inv_sqrt_3 = static_cast<value_type>(1) / std::sqrt(static_cast<value_type>(3));
1721 1000 : value_type const sqrt_2_over_sqrt_3 = std::sqrt(static_cast<value_type>(2)) / std::sqrt(static_cast<value_type>(3));
1722 1000 : r_g[0][0] = sqrt_2_over_sqrt_3;
1723 1000 : r_g[0][2] = inv_sqrt_3;
1724 1000 : r_g[2][0] = -inv_sqrt_3;
1725 1000 : r_g[2][2] = sqrt_2_over_sqrt_3;
1726 :
1727 : //std::cerr << "R_g = " << r_g << "\n";
1728 :
1729 : // $R_{rg}$ -- the product or $R_r$ and $R_g$
1730 : //
1731 2000 : matrix<T, SIZE> r_rg(r_r * r_g);
1732 :
1733 : //std::cerr << "R_rg = " << r_rg << "\n";
1734 :
1735 : // Luminance Vector
1736 : //
1737 2000 : matrix<T, SIZE> w(get_luma_vector());
1738 : //w[0][0] = RED_WEIGHT;
1739 : //w[1][0] = GREEN_WEIGHT;
1740 : //w[2][0] = BLUE_WEIGHT;
1741 : //w[3][0] = 0;
1742 :
1743 : //std::cerr << "w = " << w << "\n";
1744 :
1745 2000 : matrix<T, SIZE> l(r_rg * w);
1746 :
1747 : //std::cerr << "l = " << l << "\n";
1748 :
1749 2000 : matrix<T, SIZE> s(4, 4);
1750 1000 : s[0][2] = l[0][0] / l[2][0];
1751 1000 : s[1][2] = l[1][0] / l[2][0];
1752 :
1753 : //std::cerr << "s = " << s << "\n";
1754 :
1755 : //std::cerr << "s / s = " << (s/s) << "\n";
1756 :
1757 : //std::cerr << "M_rg * s = " << (r_rg * s) << "\n";
1758 :
1759 2000 : matrix<T, SIZE> p(r_rg);
1760 1000 : p *= s;
1761 :
1762 : // Rotate blue (rotation around Z axis)
1763 : //
1764 2000 : matrix<T, SIZE> r_b(4, 4);
1765 1000 : r_b[0][0] = rot_cos;
1766 1000 : r_b[0][1] = rot_sin;
1767 1000 : r_b[1][0] = -rot_sin;
1768 1000 : r_b[1][1] = rot_cos;
1769 :
1770 : //std::cerr << "r_b = " << r_b << "\n";
1771 :
1772 : //matrix<T, SIZE> ident(4,4);
1773 : //std::cerr
1774 : // << "---------------------------------------\n"
1775 : // << "r_rg = " << r_rg
1776 : // << "\nr_g * r_r = " << (r_g * r_r)
1777 : // << "\nr_rg^-1 = " << (ident/r_rg)
1778 : // << "\n(1/r_g/r_r) = " << (ident/r_g/r_r) << "\n"
1779 : // << "\np = " << p << "\n"
1780 : // << "\np^-1 = " << (ident/p) << "\n"
1781 : // << "---------------------------------------\n";
1782 :
1783 : // $H = R_r R_g S R_b S^{-1} R_g^{-1} R_r^{-1}$
1784 : //
1785 1000 : return *this * p * r_b / p;
1786 :
1787 : // others to test further
1788 : //return r_rg * s * r_b / s / r_rg;
1789 : //return r_rg * s * r_b / s / r_g / r_r;
1790 : //return r_r * r_g * r_b / r_g / r_r;
1791 : }
1792 : }
1793 :
1794 : #ifdef _DEBUG
1795 6000 : std::string get_last_hue_matrix()
1796 : {
1797 6000 : return f_last_hue_matrix;
1798 : }
1799 : #endif
1800 :
1801 : /** \brief Retrieve the current luma vector.
1802 : *
1803 : * By default the luma vector is set to the HDTV weights.
1804 : * This may not be desirable, though. If you would prefer to
1805 : * use NTSC (?!) or some other weights, call the
1806 : * set_luma_vector() function with the correct weights.
1807 : *
1808 : * \note
1809 : * This is often referenced as the luminance, which is not quite
1810 : * correct. What we offer here are luma weights. See:
1811 : * https://en.wikipedia.org/wiki/Luma_%28video%29 for additional
1812 : * information.
1813 : *
1814 : * You can access to the red, green, and blue weights as:
1815 : *
1816 : * \code
1817 : * matrix<T, SIZE> luma = m.get_luma_vector();
1818 : * matrix::value_type red = luma[0][0];
1819 : * matrix::value_type green = luma[1][0];
1820 : * matrix::value_type blue = luma[2][0];
1821 : * \endcode
1822 : *
1823 : * This matrix can directly be used against a 4x4 matrix.
1824 : *
1825 : * The 4th parameter is set to zero.
1826 : *
1827 : * \return A 4x1 matrix with the luma weights.
1828 : *
1829 : * \sa set_luma_vector()
1830 : */
1831 7000 : matrix<T, SIZE> get_luma_vector() const
1832 : {
1833 7000 : matrix<T, SIZE> luma(4, 1);
1834 7000 : luma[0][0] = f_luma_red;
1835 7000 : luma[1][0] = f_luma_green;
1836 7000 : luma[2][0] = f_luma_blue;
1837 7000 : return luma;
1838 : }
1839 :
1840 : /** \brief Change the luma vector.
1841 : *
1842 : * This function sets the luminance vector to the three weights
1843 : * you are passing here.
1844 : *
1845 : * The current luminance vector can be retrieved as a 4x1 matrix using
1846 : * the get_luma_vector().
1847 : *
1848 : * The recommended values are the following predefined weights:
1849 : *
1850 : * \li HDTV_RED_WEIGHT, HDTV_GREEN_WEIGHT, HDTV_BLUE_WEIGHT
1851 : * \li LED_RED_WEIGHT, LED_GREEN_WEIGHT, LED_BLUE_WEIGHT
1852 : * \li CRT_RED_WEIGHT, CRT_GREEN_WEIGHT, CRT_BLUE_WEIGHT
1853 : * \li NTSC_RED_WEIGHT, NTSC_GREEN_WEIGHT, NTSC_BLUE_WEIGHT
1854 : * \li AVERAGE_RED_WEIGHT, AVERAGE_GREEN_WEIGHT, AVERAGE_BLUE_WEIGHT
1855 : *
1856 : * The HDTV weights are those used by default if you do not call the
1857 : * set_luminance_vector() with different or even your own weights.
1858 : *
1859 : * \param[in] red_weight The luma red weight.
1860 : * \param[in] green_weight The luma green weight.
1861 : * \param[in] blue_weight The luma blue weight.
1862 : */
1863 12 : void set_luma_vector(value_type red_weight, value_type green_weight, value_type blue_weight)
1864 : {
1865 : // save the user defined weights
1866 : //
1867 12 : f_luma_red = red_weight;
1868 12 : f_luma_green = green_weight;
1869 12 : f_luma_blue = blue_weight;
1870 12 : }
1871 :
1872 : std::string to_string() const
1873 : {
1874 : std::stringstream s;
1875 :
1876 : s.precision(17);
1877 :
1878 : s << "[";
1879 : for(size_type j(0); j < f_rows; ++j)
1880 : {
1881 : s << std::endl << " [";
1882 : if(f_columns > 0)
1883 : {
1884 : s << f_vector[0 + j * f_columns];
1885 : for(size_type i(1); i < f_columns; ++i)
1886 : {
1887 : s << ", " << f_vector[i + j * f_columns];
1888 : }
1889 : }
1890 : s << "]";
1891 : }
1892 : s << std::endl << "]" << std::endl;
1893 :
1894 : return s.str();
1895 : }
1896 :
1897 : private:
1898 : friend element_ref;
1899 : friend row_ref;
1900 : friend const_row_ref;
1901 :
1902 : size_type f_rows = 0;
1903 : size_type f_columns = 0;
1904 : std::vector<T> f_vector = std::vector<T>();
1905 : value_type f_luma_red = HDTV_LUMA_RED;
1906 : value_type f_luma_green = HDTV_LUMA_GREEN;
1907 : value_type f_luma_blue = HDTV_LUMA_BLUE;
1908 : #ifdef _DEBUG
1909 : mutable std::string f_last_hue_matrix = std::string();
1910 : #endif
1911 : };
1912 :
1913 :
1914 : } // namespace snapdev
1915 :
1916 :
1917 :
1918 : /** \brief Output a matrix to a basic_ostream.
1919 : *
1920 : * This function allows one to print out a matrix. The function attempts
1921 : * to properly format the matrix in order to make it readable.
1922 : *
1923 : * \param[in] out The output stream where the matrix gets written.
1924 : * \param[in] m The actual matrix that is to be printed.
1925 : *
1926 : * \return A reference to the basic_ostream object.
1927 : */
1928 : template<class E, class S, class T, class SIZE>
1929 : std::basic_ostream<E, S> & operator << (std::basic_ostream<E, S> & out, snapdev::matrix<T, SIZE> const & m)
1930 : {
1931 : // write to a string buffer first
1932 : //
1933 : std::basic_ostringstream<E, S, std::allocator<E> > s;
1934 :
1935 : // setup the string output like the out stream
1936 : //
1937 : s.flags(out.flags());
1938 : s.imbue(out.getloc());
1939 : s.precision(out.precision());
1940 :
1941 : // output the matrix
1942 : //
1943 : s << "[";
1944 : for(typename snapdev::matrix<T, SIZE>::size_type j(0); j < m.rows(); ++j)
1945 : {
1946 : s << std::endl << " [";
1947 : if(m.columns() > 0)
1948 : {
1949 : s << static_cast<T>(m[j][0]);
1950 : for(typename snapdev::matrix<T, SIZE>::size_type i(1); i < m.columns(); ++i)
1951 : {
1952 : s << ", " << static_cast<T>(m[j][i]);
1953 : }
1954 : }
1955 : s << "]";
1956 : }
1957 : s << std::endl << "]" << std::endl;
1958 :
1959 : // buffer is ready, display in output in one go
1960 : //
1961 : return out << s.str();
1962 : }
1963 :
1964 :
1965 : // file in ve folder (matrix.c)
1966 : // http://www.graficaobscura.com/matrix/index.html
1967 : // https://ncalculators.com/matrix/3x3-matrix-multiplication-calculator.htm
1968 : // /home/alexis/m2osw/sources/freeware/libx3d/utilities/src/fast_matrix.c++
1969 :
1970 : // vim: ts=4 sw=4 et
|