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