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