LCOV - code coverage report
Current view: top level - snapdev - matrix.h (source / functions) Hit Total Coverage
Test: coverage.info Lines: 521 553 94.2 %
Date: 2021-08-22 18:14:51 Functions: 46 47 97.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13