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