PCL
Matrix.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/Matrix.h - Released 2019-11-07T10:59:34Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2019 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (http://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_Matrix_h
53 #define __PCL_Matrix_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/Container.h>
61 #include <pcl/Exception.h>
62 #include <pcl/Math.h>
63 #include <pcl/Memory.h>
64 #include <pcl/ReferenceCounter.h>
65 #include <pcl/Rotate.h> // pcl::Reverse()
66 #include <pcl/Utility.h>
67 #include <pcl/Vector.h>
68 
69 #ifndef __PCL_NO_MATRIX_STATISTICS
70 # include <pcl/Selection.h>
71 # include <pcl/Sort.h>
72 #endif
73 
74 #if !defined( __PCL_NO_MATRIX_IMAGE_RENDERING ) && !defined( __PCL_NO_MATRIX_IMAGE_CONVERSION )
75 # include <pcl/Image.h>
76 # include <pcl/ImageVariant.h>
77 #endif
78 
79 #if !defined( __PCL_NO_MATRIX_PHASE_MATRICES ) && !defined( __PCL_NO_VECTOR_INSTANTIATE )
80 # include <pcl/Complex.h>
81 #endif
82 
83 /*
84  * Valid filter kernel sizes are odd integers >= 3. This macro is used for the
85  * KernelFilter and SeparableFilter classes to ensure validity of filter
86  * matrices and vectors.
87  */
88 #define PCL_VALID_KERNEL_SIZE( n ) (n ? Max( 3, n|1 ) : 0)
89 
90 namespace pcl
91 {
92 
93 // ----------------------------------------------------------------------------
94 
121 template <typename T>
122 class PCL_CLASS GenericMatrix : public DirectContainer<T>
123 {
124 public:
125 
129  typedef T element;
130 
135 
141  typedef element* block_iterator;
142 
148  typedef const element* const_block_iterator;
149 
155  {
156  m_data = new Data;
157  }
158 
168  GenericMatrix( int rows, int cols )
169  {
170  m_data = new Data( rows, cols );
171  }
172 
180  GenericMatrix( const element& x, int rows, int cols )
181  {
182  m_data = new Data( rows, cols );
183  pcl::Fill( m_data->Begin(), m_data->End(), x );
184  }
185 
199  template <typename T1>
200  GenericMatrix( const T1* a, int rows, int cols )
201  {
202  m_data = new Data( rows, cols );
203  if ( a != nullptr )
204  {
205  block_iterator c = m_data->Begin();
206  for ( const T1* b = a + NumberOfElements(); a < b; ++a, ++c )
207  *c = element( *a );
208  }
209  }
210 
222  template <typename T1>
223  GenericMatrix( const T1& a00, const T1& a01, const T1& a02,
224  const T1& a10, const T1& a11, const T1& a12,
225  const T1& a20, const T1& a21, const T1& a22 )
226  {
227  m_data = new Data( 3, 3 );
228  block_iterator v = m_data->Begin();
229  v[0] = element( a00 ); v[1] = element( a01 ); v[2] = element( a02 );
230  v[3] = element( a10 ); v[4] = element( a11 ); v[5] = element( a12 );
231  v[6] = element( a20 ); v[7] = element( a21 ); v[8] = element( a22 );
232  }
233 
238  GenericMatrix( const GenericMatrix& x ) : m_data( x.m_data )
239  {
240  m_data->Attach();
241  }
242 
246  GenericMatrix( GenericMatrix&& x ) : m_data( x.m_data )
247  {
248  x.m_data = nullptr;
249  }
250 
273  GenericMatrix( const GenericMatrix& x, int i0, int j0, int rows, int cols )
274  {
275  i0 = Range( i0, 0, Max( 0, x.Rows()-1 ) );
276  j0 = Range( j0, 0, Max( 0, x.Cols()-1 ) );
277  rows = Range( rows, 0, Max( 0, x.Rows()-i0 ) );
278  cols = Range( cols, 0, Max( 0, x.Cols()-j0 ) );
279  m_data = new Data( rows, cols );
280  for ( int i = 0; i < m_data->Rows(); ++i, ++i0 )
281  for ( int j = 0; j < m_data->Cols(); ++j, ++j0 )
282  m_data->v[i][j] = x.m_data->v[i0][j0];
283  }
284 
285 #ifndef __PCL_NO_MATRIX_IMAGE_CONVERSION
286 
323  template <class P>
324  GenericMatrix( const GenericImage<P>& image, const Rect& rect = Rect( 0 ), int channel = -1 )
325  {
326  GenericMatrix M = FromImage( image, rect, channel );
327  pcl::Swap( m_data, M.m_data );
328  }
329 
344  GenericMatrix( const ImageVariant& image, const Rect& rect = Rect( 0 ), int channel = -1 )
345  {
346  GenericMatrix M = FromImage( image, rect, channel );
347  pcl::Swap( m_data, M.m_data );
348  }
349 
350 #endif // !__PCL_NO_MATRIX_IMAGE_CONVERSION
351 
357  virtual ~GenericMatrix()
358  {
359  if ( m_data != nullptr )
360  {
361  DetachFromData();
362  m_data = nullptr;
363  }
364  }
365 
369  void Clear()
370  {
371  if ( !IsEmpty() )
372  if ( m_data->IsUnique() )
373  m_data->Deallocate();
374  else
375  {
376  Data* newData = new Data( 0, 0 );
377  DetachFromData();
378  m_data = newData;
379  }
380  }
381 
394  GenericMatrix& operator =( const GenericMatrix& x )
395  {
396  Assign( x );
397  return *this;
398  }
399 
412  void Assign( const GenericMatrix& x )
413  {
414  x.m_data->Attach();
415  DetachFromData();
416  m_data = x.m_data;
417  }
418 
422  GenericMatrix& operator =( GenericMatrix&& x )
423  {
424  Transfer( x );
425  return *this;
426  }
427 
441  {
442  DetachFromData();
443  m_data = x.m_data;
444  x.m_data = nullptr;
445  }
446 
460  {
461  DetachFromData();
462  m_data = x.m_data;
463  x.m_data = nullptr;
464  }
465 
472  friend void Swap( GenericMatrix& x1, GenericMatrix& x2 )
473  {
474  pcl::Swap( x1.m_data, x2.m_data );
475  }
476 
485  GenericMatrix& operator =( const element& x )
486  {
487  if ( !IsUnique() )
488  {
489  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
490  DetachFromData();
491  m_data = newData;
492  }
493  pcl::Fill( m_data->Begin(), m_data->End(), x );
494  return *this;
495  }
496 
497 #define IMPLEMENT_SCALAR_ASSIGN_OP( op ) \
498  if ( IsUnique() ) \
499  { \
500  block_iterator a = m_data->Begin(); \
501  const_block_iterator b = m_data->End(); \
502  for ( ; a < b; ++a ) \
503  *a op##= x; \
504  } \
505  else \
506  { \
507  Data* newData = new Data( m_data->Rows(), m_data->Cols() ); \
508  block_iterator a = newData->Begin(); \
509  const_block_iterator b = newData->End(); \
510  const_block_iterator c = m_data->Begin(); \
511  for ( ; a < b; ++a, ++c ) \
512  *a = *c op x; \
513  DetachFromData(); \
514  m_data = newData; \
515  } \
516  return *this;
517 
526  GenericMatrix& operator +=( const element& x )
527  {
528  IMPLEMENT_SCALAR_ASSIGN_OP( + )
529  }
530 
539  GenericMatrix& operator -=( const element& x )
540  {
541  IMPLEMENT_SCALAR_ASSIGN_OP( - )
542  }
543 
552  GenericMatrix& operator *=( const element& x )
553  {
554  IMPLEMENT_SCALAR_ASSIGN_OP( * )
555  }
556 
565  GenericMatrix& operator /=( const element& x )
566  {
567  IMPLEMENT_SCALAR_ASSIGN_OP( / )
568  }
569 
578  GenericMatrix& operator ^=( const element& x )
579  {
580  if ( IsUnique() )
581  {
582  block_iterator a = m_data->Begin();
583  const_block_iterator b = m_data->End();
584  for ( ; a < b; ++a )
585  *a = pcl::Pow( *a, x );
586  }
587  else
588  {
589  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
590  block_iterator a = newData->Begin();
591  const_block_iterator b = newData->End();
592  const_block_iterator c = m_data->Begin();
593  for ( ; a < b; ++a, ++c )
594  *a = pcl::Pow( *c, x );
595  DetachFromData();
596  m_data = newData;
597  }
598  return *this;
599  }
600 
601 #undef IMPLEMENT_SCALAR_ASSIGN_OP
602 
603 #ifndef __PCL_NO_MATRIX_ELEMENT_WISE_ARITHMETIC_OPERATIONS
604 
605 #define IMPLEMENT_ELEMENT_WISE_ASSIGN_OP( op ) \
606  if ( IsUnique() ) \
607  { \
608  block_iterator a = m_data->Begin(); \
609  const_block_iterator b = pcl::Min( m_data->End(), a + x.NumberOfElements() ); \
610  const_block_iterator c = x.m_data->Begin(); \
611  for ( ; a < b; ++a, ++c ) \
612  *a op##= *c; \
613  } \
614  else \
615  { \
616  Data* newData = new Data( m_data->Rows(), m_data->Cols() ); \
617  block_iterator a = newData->Begin(); \
618  const_block_iterator b = pcl::Min( newData->End(), a + x.NumberOfElements() ); \
619  const_block_iterator c = m_data->Begin(); \
620  const_block_iterator d = x.m_data->Begin(); \
621  for ( ; a < b; ++a, ++c, ++d ) \
622  *a = *c op *d; \
623  DetachFromData(); \
624  m_data = newData; \
625  } \
626  return *this;
627 
640  GenericMatrix& operator +=( const GenericMatrix& x )
641  {
642  IMPLEMENT_ELEMENT_WISE_ASSIGN_OP( + )
643  }
644 
657  GenericMatrix& operator -=( const GenericMatrix& x )
658  {
659  IMPLEMENT_ELEMENT_WISE_ASSIGN_OP( - )
660  }
661 
674  GenericMatrix& operator *=( const GenericMatrix& x )
675  {
676  IMPLEMENT_ELEMENT_WISE_ASSIGN_OP( * )
677  }
678 
691  GenericMatrix& operator /=( const GenericMatrix& x )
692  {
693  IMPLEMENT_ELEMENT_WISE_ASSIGN_OP( / )
694  }
695 
707  GenericMatrix& operator ^=( const GenericMatrix& x )
708  {
709  if ( IsUnique() )
710  {
711  block_iterator a = m_data->Begin();
712  const_block_iterator b = pcl::Min( m_data->End(), a + x.NumberOfElements() );
713  const_block_iterator c = x.m_data->Begin();
714  for ( ; a < b; ++a, ++c )
715  *a = pcl::Pow( *a, *c );
716  }
717  else
718  {
719  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
720  block_iterator a = newData->Begin();
721  const_block_iterator b = pcl::Min( newData->End(), a + x.NumberOfElements() );
722  const_block_iterator c = m_data->Begin();
723  const_block_iterator d = x.m_data->Begin();
724  for ( ; a < b; ++a, ++c, ++d )
725  *a = pcl::Pow( *c, *d );
726  DetachFromData();
727  m_data = newData;
728  }
729  return *this;
730  }
731 
732 #undef IMPLEMENT_ELEMENT_WISE_ASSIGN_OP
733 
734 #endif // __PCL_NO_MATRIX_ELEMENT_WISE_ARITHMETIC_OPERATIONS
735 
742  {
743  GenericMatrix R( m_data->Rows(), m_data->Cols() );
744  block_iterator r = R.m_data->Begin();
745  const_block_iterator s = R.m_data->End();
746  const_block_iterator a = m_data->Begin();
747  for ( ; r < s; ++r, ++a )
748  *r = *a * *a;
749  return R;
750  }
751 
758  void SetSqr()
759  {
760  if ( IsUnique() )
761  {
762  block_iterator a = m_data->Begin();
763  const_block_iterator b = m_data->End();
764  for ( ; a < b; ++a )
765  *a *= *a;
766  }
767  else
768  {
769  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
770  block_iterator a = newData->Begin();
771  const_block_iterator b = newData->End();
772  const_block_iterator c = m_data->Begin();
773  for ( ; a < b; ++a, ++c )
774  *a = *c * *c;
775  DetachFromData();
776  m_data = newData;
777  }
778  }
779 
786  {
787  GenericMatrix R( m_data->Rows(), m_data->Cols() );
788  block_iterator r = R.m_data->Begin();
789  const_block_iterator s = R.m_data->End();
790  const_block_iterator a = m_data->Begin();
791  for ( ; r < s; ++r, ++a )
792  *r = pcl::Sqrt( *a );
793  return R;
794  }
795 
802  void SetSqrt()
803  {
804  if ( IsUnique() )
805  {
806  block_iterator a = m_data->Begin();
807  const_block_iterator b = m_data->End();
808  for ( ; a < b; ++a )
809  *a = pcl::Sqrt( *a );
810  }
811  else
812  {
813  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
814  block_iterator a = newData->Begin();
815  const_block_iterator b = newData->End();
816  const_block_iterator c = m_data->Begin();
817  for ( ; a < b; ++a, ++c )
818  *a = pcl::Sqrt( *c );
819  DetachFromData();
820  m_data = newData;
821  }
822  }
823 
830  {
831  GenericMatrix R( m_data->Rows(), m_data->Cols() );
832  block_iterator r = R.m_data->Begin();
833  const_block_iterator s = R.m_data->End();
834  const_block_iterator a = m_data->Begin();
835  for ( ; r < s; ++r, ++a )
836  *r = pcl::Abs( *a );
837  return R;
838  }
839 
846  void SetAbs()
847  {
848  if ( IsUnique() )
849  {
850  block_iterator a = m_data->Begin();
851  const_block_iterator b = m_data->End();
852  for ( ; a < b; ++a )
853  *a = pcl::Abs( *a );
854  }
855  else
856  {
857  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
858  block_iterator a = newData->Begin();
859  const_block_iterator b = newData->End();
860  const_block_iterator c = m_data->Begin();
861  for ( ; a < b; ++a, ++c )
862  *a = pcl::Abs( *c );
863  DetachFromData();
864  m_data = newData;
865  }
866  }
867 
871  bool IsUnique() const
872  {
873  return m_data->IsUnique();
874  }
875 
880  bool IsAliasOf( const GenericMatrix& x ) const
881  {
882  return m_data == x.m_data;
883  }
884 
893  {
894  if ( !IsUnique() )
895  {
896  Data* newData = new Data( m_data->Rows(), m_data->Cols() );
897  pcl::Copy( newData->Begin(), m_data->Begin(), m_data->End() );
898  DetachFromData();
899  m_data = newData;
900  }
901  }
902 
907  int Rows() const
908  {
909  return m_data->Rows();
910  }
911 
916  int Cols() const
917  {
918  return m_data->Cols();
919  }
920 
927  int Columns() const
928  {
929  return Cols();
930  }
931 
947  bool IsValid() const
948  {
949  return m_data != nullptr;
950  }
951 
956  bool IsEmpty() const
957  {
958  return Rows() == 0 || Cols() == 0;
959  }
960 
966  operator bool() const
967  {
968  return !IsEmpty();
969  }
970 
976  {
977  return m_data->NumberOfElements();
978  }
979 
984  size_type Size() const
985  {
986  return m_data->Size();
987  }
988 
996  bool operator ==( const GenericMatrix& x ) const
997  {
998  return IsAliasOf( x ) || SameDimensions( x ) && pcl::Equal( Begin(), x.Begin(), x.End() );
999  }
1000 
1010  bool operator <( const GenericMatrix& x ) const
1011  {
1012  return !IsAliasOf( x ) && pcl::Compare( Begin(), End(), x.Begin(), x.End() ) < 0;
1013  }
1014 
1019  bool SameDimensions( const GenericMatrix& x ) const
1020  {
1021  return Rows() == x.Rows() && Cols() == x.Cols();
1022  }
1023 
1033  bool SameElements( const GenericMatrix& x ) const
1034  {
1035  if ( IsAliasOf( x ) )
1036  return true;
1037  if ( NumberOfElements() != x.NumberOfElements() )
1038  return false;
1039  const_block_iterator a = Begin();
1040  const_block_iterator b = End();
1041  const_block_iterator c = x.Begin();
1042  for ( ; a < b; ++a, ++c )
1043  if ( *a != *c )
1044  return false;
1045  return true;
1046  }
1047 
1054  element& Element( int i, int j )
1055  {
1056  EnsureUnique();
1057  return m_data->Element( i, j );
1058  }
1059 
1064  const element& Element( int i, int j ) const
1065  {
1066  return m_data->Element( i, j );
1067  }
1068 
1078  block_iterator operator []( int i )
1079  {
1080  EnsureUnique();
1081  return m_data->v[i];
1082  }
1083 
1090  const_block_iterator operator []( int i ) const
1091  {
1092  return m_data->v[i];
1093  }
1094 
1107  block_iterator Begin()
1108  {
1109  EnsureUnique();
1110  return m_data->Begin();
1111  }
1112 
1122  const_block_iterator Begin() const
1123  {
1124  return m_data->Begin();
1125  }
1126 
1130  const_block_iterator ConstBegin() const
1131  {
1132  return Begin();
1133  }
1134 
1141  block_iterator operator *()
1142  {
1143  return Begin();
1144  }
1145 
1152  const_block_iterator operator *() const
1153  {
1154  return Begin();
1155  }
1156 
1169  block_iterator End()
1170  {
1171  EnsureUnique();
1172  return m_data->End();
1173  }
1174 
1185  const_block_iterator End() const
1186  {
1187  return m_data->End();
1188  }
1189 
1193  const_block_iterator ConstEnd() const
1194  {
1195  return End();
1196  }
1197 
1216  block_iterator* DataPtr()
1217  {
1218  return m_data->v;
1219  }
1220 
1231  block_iterator RowPtr( int i )
1232  {
1233  return m_data->v[i];
1234  }
1235 
1236 #ifndef __PCL_NO_STL_COMPATIBLE_ITERATORS
1237 
1240  block_iterator begin()
1241  {
1242  return Begin();
1243  }
1244 
1248  const_block_iterator begin() const
1249  {
1250  return Begin();
1251  }
1252 
1256  block_iterator end()
1257  {
1258  return End();
1259  }
1260 
1264  const_block_iterator end() const
1265  {
1266  return End();
1267  }
1268 #endif // !__PCL_NO_STL_COMPATIBLE_ITERATORS
1269 
1273  vector RowVector( int i ) const
1274  {
1275  vector r( m_data->Cols() );
1276  for ( int j = 0; j < m_data->Cols(); ++j )
1277  r[j] = m_data->v[i][j];
1278  return r;
1279  }
1280 
1284  vector ColumnVector( int j ) const
1285  {
1286  vector c( m_data->Rows() );
1287  for ( int i = 0; i < m_data->Rows(); ++i )
1288  c[i] = m_data->v[i][j];
1289  return c;
1290  }
1291 
1297  vector ColVector( int j ) const
1298  {
1299  return ColumnVector( j );
1300  }
1301 
1306  template <class V>
1307  void SetRow( int i, const V& r )
1308  {
1309  EnsureUnique();
1310  for ( int j = 0; j < m_data->Cols() && j < r.Length(); ++j )
1311  m_data->v[i][j] = element( r[j] );
1312  }
1313 
1318  template <class V>
1319  void SetColumn( int j, const V& c )
1320  {
1321  EnsureUnique();
1322  for ( int i = 0; i < m_data->Rows() && i < c.Length(); ++i )
1323  m_data->v[i][j] = element( c[i] );
1324  }
1325 
1332  template <class V>
1333  void SetCol( int j, const V& c )
1334  {
1335  SetColumn( j, c );
1336  }
1337 
1343  static GenericMatrix FromRowVector( const vector& r )
1344  {
1345  GenericMatrix R( element( 0 ), r.Length(), r.Length() );
1346  for ( int j = 0; j < r.Length(); ++j )
1347  R.m_data->v[0][j] = r[j];
1348  return R;
1349  }
1350 
1356  static GenericMatrix FromColumnVector( const vector& c )
1357  {
1358  GenericMatrix C( element( 0 ), c.Length(), c.Length() );
1359  for ( int i = 0; i < c.Length(); ++i )
1360  C.m_data->v[i][0] = c[i];
1361  return C;
1362  }
1363 
1367  static GenericMatrix FromColVector( const vector& c )
1368  {
1369  return FromColumnVector( c );
1370  }
1371 
1378  static GenericMatrix UnitMatrix( int n )
1379  {
1380  GenericMatrix One( element( 0 ), n, n );
1381  for ( int i = 0; i < n; ++i )
1382  One[i][i] = element( 1 );
1383  return One;
1384  }
1385 
1393  {
1394  GenericMatrix Tr( m_data->Cols(), m_data->Rows() );
1395  for ( int i = 0; i < m_data->Rows(); ++i )
1396  for ( int j = 0; j < m_data->Cols(); ++j )
1397  Tr.m_data->v[j][i] = m_data->v[i][j];
1398  return Tr;
1399  }
1400 
1411  GenericMatrix Inverse() const;
1412 
1423  void Invert();
1424 
1432  void Flip()
1433  {
1434  EnsureUnique();
1435  pcl::Reverse( m_data->Begin(), m_data->End() );
1436  }
1437 
1443  {
1444  GenericMatrix Tf( m_data->Cols(), m_data->Rows() );
1445  pcl::CopyReversed( Tf.m_data->End(), m_data->Begin(), m_data->End() );
1446  return Tf;
1447  }
1448 
1472  void RotateX( double sphi, double cphi )
1473  {
1474  PCL_PRECONDITION( Rows() == 3 && Cols() == 3 )
1475  block_iterator A1 = m_data->v[1];
1476  block_iterator A2 = m_data->v[2];
1477  double a10 = cphi*A1[0] + sphi*A2[0];
1478  double a11 = cphi*A1[1] + sphi*A2[1];
1479  double a12 = cphi*A1[2] + sphi*A2[2];
1480  double a20 = -sphi*A1[0] + cphi*A2[0];
1481  double a21 = -sphi*A1[1] + cphi*A2[1];
1482  double a22 = -sphi*A1[2] + cphi*A2[2];
1483  A1[0] = element( a10 );
1484  A1[1] = element( a11 );
1485  A1[2] = element( a12 );
1486  A2[0] = element( a20 );
1487  A2[1] = element( a21 );
1488  A2[2] = element( a22 );
1489  }
1490 
1499  void RotateX( double phi )
1500  {
1501  double sphi, cphi;
1502  SinCos( phi, sphi, cphi );
1503  RotateX( sphi, cphi );
1504  }
1505 
1511  GenericMatrix RotatedX( double sphi, double cphi ) const
1512  {
1513  GenericMatrix R( *this );
1514  R.RotateX( sphi, cphi );
1515  return R;
1516  }
1517 
1522  GenericMatrix RotatedX( double phi ) const
1523  {
1524  GenericMatrix R( *this );
1525  R.RotateX( phi );
1526  return R;
1527  }
1528 
1552  void RotateY( double sphi, double cphi )
1553  {
1554  PCL_PRECONDITION( Rows() == 3 && Cols() == 3 )
1555  block_iterator A0 = m_data->v[0];
1556  block_iterator A2 = m_data->v[2];
1557  double a00 = cphi*A0[0] - sphi*A2[0];
1558  double a01 = cphi*A0[1] - sphi*A2[1];
1559  double a02 = cphi*A0[2] - sphi*A2[2];
1560  double a20 = sphi*A0[0] + cphi*A2[0];
1561  double a21 = sphi*A0[1] + cphi*A2[1];
1562  double a22 = sphi*A0[2] + cphi*A2[2];
1563  A0[0] = a00;
1564  A0[1] = a01;
1565  A0[2] = a02;
1566  A2[0] = a20;
1567  A2[1] = a21;
1568  A2[2] = a22;
1569  }
1570 
1579  void RotateY( double phi )
1580  {
1581  double sphi, cphi;
1582  SinCos( phi, sphi, cphi );
1583  RotateY( sphi, cphi );
1584  }
1585 
1591  GenericMatrix RotatedY( double sphi, double cphi ) const
1592  {
1593  GenericMatrix R( *this );
1594  R.RotateY( sphi, cphi );
1595  return R;
1596  }
1597 
1602  GenericMatrix RotatedY( double phi ) const
1603  {
1604  GenericMatrix R( *this );
1605  R.RotateY( phi );
1606  return R;
1607  }
1608 
1632  void RotateZ( double sphi, double cphi )
1633  {
1634  PCL_PRECONDITION( Rows() == 3 && Cols() == 3 )
1635  block_iterator A0 = m_data->v[0];
1636  block_iterator A1 = m_data->v[1];
1637  double a00 = cphi*A0[0] + sphi*A1[0];
1638  double a01 = cphi*A0[1] + sphi*A1[1];
1639  double a02 = cphi*A0[2] + sphi*A1[2];
1640  double a10 = -sphi*A0[0] + cphi*A1[0];
1641  double a11 = -sphi*A0[1] + cphi*A1[1];
1642  double a12 = -sphi*A0[2] + cphi*A1[2];
1643  A0[0] = a00;
1644  A0[1] = a01;
1645  A0[2] = a02;
1646  A1[0] = a10;
1647  A1[1] = a11;
1648  A1[2] = a12;
1649  }
1650 
1659  void RotateZ( double phi )
1660  {
1661  double sphi, cphi;
1662  SinCos( phi, sphi, cphi );
1663  RotateZ( sphi, cphi );
1664  }
1665 
1671  GenericMatrix RotatedZ( double sphi, double cphi ) const
1672  {
1673  GenericMatrix R( *this );
1674  R.RotateZ( sphi, cphi );
1675  return R;
1676  }
1677 
1682  GenericMatrix RotatedZ( double phi ) const
1683  {
1684  GenericMatrix R( *this );
1685  R.RotateZ( phi );
1686  return R;
1687  }
1688 
1701  void Truncate( const element& f0 = element( 0 ), const element& f1 = element( 1 ) )
1702  {
1703  EnsureUnique();
1704  block_iterator a = m_data->Begin();
1705  block_iterator b = m_data->End();
1706  for ( ; a < b; ++a )
1707  if ( *a < f0 )
1708  *a = f0;
1709  else if ( *a > f1 )
1710  *a = f1;
1711  }
1712 
1716  GenericMatrix Truncated( const element& f0 = element( 0 ), const element& f1 = element( 1 ) ) const
1717  {
1718  GenericMatrix R( *this );
1719  R.Truncate( f0, f1 );
1720  return R;
1721  }
1722 
1738  void Rescale( const element& f0 = element( 0 ), const element& f1 = element( 1 ) )
1739  {
1740  element v0 = MinElement();
1741  element v1 = MaxElement();
1742  if ( v0 != f0 || v1 != f1 )
1743  {
1744  EnsureUnique();
1745  if ( v0 != v1 )
1746  {
1747  if ( f0 != f1 )
1748  {
1749  block_iterator a = m_data->Begin();
1750  block_iterator b = m_data->End();
1751  double d = (double( f1 ) - double( f0 ))/(double( v1 ) - double( v0 ));
1752  if ( f0 == element( 0 ) )
1753  for ( ; a < b; ++a )
1754  *a = element( d*(*a - v0) );
1755  else
1756  for ( ; a < b; ++a )
1757  *a = element( d*(*a - v0) + f0 );
1758  }
1759  else
1760  pcl::Fill( m_data->Begin(), m_data->End(), f0 );
1761  }
1762  else
1763  pcl::Fill( m_data->Begin(), m_data->End(), pcl::Range( v0, f0, f1 ) );
1764  }
1765  }
1766 
1770  GenericMatrix Rescaled( const element& f0 = element( 0 ), const element& f1 = element( 1 ) ) const
1771  {
1772  GenericMatrix R( *this );
1773  R.Rescale( f0, f1 );
1774  return R;
1775  }
1776 
1780  void Sort()
1781  {
1782  EnsureUnique();
1783  pcl::Sort( m_data->Begin(), m_data->End() );
1784  }
1785 
1790  {
1791  GenericMatrix R( *this );
1792  R.Sort();
1793  return R;
1794  }
1795 
1800  {
1801  EnsureUnique();
1802  pcl::Sort( m_data->Begin(), m_data->End(),
1803  []( const element& a, const element& b ){ return b < a; } );
1804  }
1805 
1810  {
1811  GenericMatrix R( *this );
1812  R.ReverseSort();
1813  return R;
1814  }
1815 
1821  template <class BP>
1822  void Sort( BP p )
1823  {
1824  EnsureUnique();
1825  pcl::Sort( m_data->Begin(), m_data->End(), p );
1826  }
1827 
1832  template <class BP>
1833  GenericMatrix Sorted( BP p ) const
1834  {
1835  GenericMatrix R( *this );
1836  R.Sort( p );
1837  return R;
1838  }
1839 
1844  const_block_iterator Find( const element& x ) const
1845  {
1846  const_block_iterator p = pcl::LinearSearch( m_data->Begin(), m_data->End(), x );
1847  return (p != m_data->End()) ? p : 0;
1848  }
1849 
1855  const_block_iterator FindFirst( const element& x ) const
1856  {
1857  return Find( x );
1858  }
1859 
1864  const_block_iterator FindLast( const element& x ) const
1865  {
1866  const_block_iterator p = pcl::LinearSearchLast( m_data->Begin(), m_data->End(), x );
1867  return (p != m_data->End()) ? p : -1;
1868  }
1869 
1873  bool Contains( const element& x ) const
1874  {
1875  return pcl::LinearSearch( m_data->Begin(), m_data->End(), x ) != m_data->End();
1876  }
1877 
1878 #ifndef __PCL_NO_MATRIX_STATISTICS
1879 
1884  element MinElement() const
1885  {
1886  if ( !IsEmpty() )
1887  return *MinItem( m_data->Begin(), m_data->End() );
1888  return element( 0 );
1889  }
1890 
1897  element MinElement( Point& p ) const
1898  {
1899  if ( !IsEmpty() )
1900  {
1901  const_block_iterator m = MinItem( m_data->Begin(), m_data->End() );
1902  distance_type d = m - m_data->Begin();
1903  p.y = int( d/m_data->Cols() );
1904  p.x = int( d%m_data->Cols() );
1905  return *m;
1906  }
1907  p = 0;
1908  return element( 0 );
1909  }
1910 
1915  element MaxElement() const
1916  {
1917  if ( !IsEmpty() )
1918  return *MaxItem( m_data->Begin(), m_data->End() );
1919  return element( 0 );
1920  }
1921 
1928  element MaxElement( Point& p ) const
1929  {
1930  if ( !IsEmpty() )
1931  {
1932  const_block_iterator m = MaxItem( m_data->Begin(), m_data->End() );
1933  int d = m - m_data->Begin();
1934  p.y = d/m_data->Cols();
1935  p.x = d%m_data->Cols();
1936  return *m;
1937  }
1938  p = 0;
1939  return element( 0 );
1940  }
1941 
1947  double Sum() const
1948  {
1949  return pcl::Sum( m_data->Begin(), m_data->End() );
1950  }
1951 
1958  double StableSum() const
1959  {
1960  return pcl::StableSum( m_data->Begin(), m_data->End() );
1961  }
1962 
1968  double Modulus() const
1969  {
1970  return pcl::Modulus( m_data->Begin(), m_data->End() );
1971  }
1972 
1979  double StableModulus() const
1980  {
1981  return pcl::StableModulus( m_data->Begin(), m_data->End() );
1982  }
1983 
1989  double SumOfSquares() const
1990  {
1991  return pcl::SumOfSquares( m_data->Begin(), m_data->End() );
1992  }
1993 
2000  double StableSumOfSquares() const
2001  {
2002  return pcl::StableSumOfSquares( m_data->Begin(), m_data->End() );
2003  }
2004 
2010  double Mean() const
2011  {
2012  return pcl::Mean( m_data->Begin(), m_data->End() );
2013  }
2014 
2021  double StableMean() const
2022  {
2023  return pcl::StableMean( m_data->Begin(), m_data->End() );
2024  }
2025 
2031  double Variance() const
2032  {
2033  return pcl::Variance( m_data->Begin(), m_data->End() );
2034  }
2035 
2042  double StdDev() const
2043  {
2044  return pcl::StdDev( m_data->Begin(), m_data->End() );
2045  }
2046 
2059  double Median()
2060  {
2061  EnsureUnique();
2062  return pcl::Median( m_data->Begin(), m_data->End() );
2063  }
2064 
2077  double Median() const
2078  {
2079  return GenericMatrix( *this ).Median();
2080  }
2081 
2096  double AvgDev( double center ) const
2097  {
2098  return pcl::AvgDev( m_data->Begin(), m_data->End(), center );
2099  }
2100 
2116  double StableAvgDev( double center ) const
2117  {
2118  return pcl::StableAvgDev( m_data->Begin(), m_data->End(), center );
2119  }
2120 
2133  double AvgDev() const
2134  {
2135  return pcl::AvgDev( m_data->Begin(), m_data->End() );
2136  }
2137 
2151  double StableAvgDev() const
2152  {
2153  return pcl::StableAvgDev( m_data->Begin(), m_data->End() );
2154  }
2155 
2167  double MAD( double center ) const
2168  {
2169  return pcl::MAD( m_data->Begin(), m_data->End(), center );
2170  }
2171 
2182  double MAD() const
2183  {
2184  return pcl::MAD( m_data->Begin(), m_data->End() );
2185  }
2186 
2211  double BiweightMidvariance( double center, double sigma, int k = 9 ) const
2212  {
2213  return pcl::BiweightMidvariance( m_data->Begin(), m_data->End(), center, sigma, k );
2214  }
2215 
2234  double BiweightMidvariance( int k = 9 ) const
2235  {
2236  double center = Median();
2237  double sigma = 1.4826*MAD( center );
2238  return pcl::BiweightMidvariance( m_data->Begin(), m_data->End(), center, sigma, k );
2239  }
2240 
2263  double BendMidvariance( double center, double beta = 0.2 ) const
2264  {
2265  return pcl::BendMidvariance( m_data->Begin(), m_data->End(), center, beta );
2266  }
2267 
2287  double BendMidvariance( double beta = 0.2 ) const
2288  {
2289  return pcl::BendMidvariance( m_data->Begin(), m_data->End(), Median(), beta );
2290  }
2291 
2313  double Sn() const
2314  {
2315  GenericMatrix A( *this );
2316  return pcl::Sn( A.Begin(), A.End() );
2317  }
2318 
2339  double Qn() const
2340  {
2341  GenericMatrix A( *this );
2342  return pcl::Qn( A.Begin(), A.End() );
2343  }
2344 
2345 #endif // !__PCL_NO_MATRIX_STATISTICS
2346 
2355  uint64 Hash64( uint64 seed = 0 ) const
2356  {
2357  return pcl::Hash64( m_data->Begin(), m_data->Size(), seed );
2358  }
2359 
2368  uint32 Hash32( uint32 seed = 0 ) const
2369  {
2370  return pcl::Hash32( m_data->Begin(), m_data->Size(), seed );
2371  }
2372 
2377  uint64 Hash( uint64 seed = 0 ) const
2378  {
2379  return Hash64( seed );
2380  }
2381 
2398  template <class S, typename SP>
2399  S& ToSeparated( S& s, SP separator ) const
2400  {
2401  const_block_iterator i = m_data->Begin(), j = m_data->End();
2402  if ( i < j )
2403  {
2404  s.Append( S( *i ) );
2405  if ( ++i < j )
2406  do
2407  {
2408  s.Append( separator );
2409  s.Append( S( *i ) );
2410  }
2411  while ( ++i < j );
2412  }
2413  return s;
2414  }
2415 
2438  template <class S, typename SP, class AF>
2439  S& ToSeparated( S& s, SP separator, AF append ) const
2440  {
2441  const_block_iterator i = m_data->Begin(), j = m_data->End();
2442  if ( i < j )
2443  {
2444  append( s, S( *i ) );
2445  if ( ++i < j )
2446  {
2447  S p( separator );
2448  do
2449  {
2450  append( s, p );
2451  append( s, S( *i ) );
2452  }
2453  while ( ++i < j );
2454  }
2455  }
2456  return s;
2457  }
2458 
2467  template <class S>
2468  S& ToCommaSeparated( S& s ) const
2469  {
2470  return ToSeparated( s, ',' );
2471  }
2472 
2481  template <class S>
2482  S& ToSpaceSeparated( S& s ) const
2483  {
2484  return ToSeparated( s, ' ' );
2485  }
2486 
2495  template <class S>
2496  S& ToTabSeparated( S& s ) const
2497  {
2498  return ToSeparated( s, '\t' );
2499  }
2500 
2501 #ifndef __PCL_NO_MATRIX_PHASE_MATRICES
2502 
2508  {
2509  if ( B.Rows() != A.Rows() || B.Cols() != A.Cols() )
2510  throw Error( "Invalid matrix dimensions in PhaseCorrelationMatrix()" );
2511  GenericMatrix R( A.Rows(), A.Cols() );
2512  PhaseCorrelationMatrix( R.Begin(), R.End(), A.Begin(), B.Begin() );
2513  return R;
2514  }
2515 
2521  {
2522  if ( B.Rows() != A.Rows() || B.Cols() != A.Cols() )
2523  throw Error( "Invalid matrix dimensions in CrossPowerSpectrumMatrix()" );
2524  GenericMatrix R( A.Rows(), A.Cols() );
2525  CrossPowerSpectrumMatrix( R.Begin(), R.End(), A.Begin(), B.Begin() );
2526  return R;
2527  }
2528 
2529 #endif // !__PCL_NO_MATRIX_PHASE_MATRICES
2530 
2531 #ifndef __PCL_NO_MATRIX_IMAGE_RENDERING
2532 
2551  template <class P>
2552  void ToImage( GenericImage<P>& image ) const
2553  {
2554  image.AllocateData( Cols(), Rows() );
2555  typename P::sample* v = *image;
2556  const_block_iterator a = m_data->Begin();
2557  const_block_iterator b = m_data->End();
2558  for ( ; a < b; ++a, ++v )
2559  *v = P::ToSample( *a );
2560  }
2561 
2577  void ToImage( ImageVariant& image ) const
2578  {
2579  if ( !image )
2580  image.CreateFloatImage();
2581 
2582  if ( image.IsFloatSample() )
2583  switch ( image.BitsPerSample() )
2584  {
2585  case 32: ToImage( static_cast<Image&>( *image ) ); break;
2586  case 64: ToImage( static_cast<DImage&>( *image ) ); break;
2587  }
2588  else if ( image.IsComplexSample() )
2589  switch ( image.BitsPerSample() )
2590  {
2591  case 32: ToImage( static_cast<ComplexImage&>( *image ) ); break;
2592  case 64: ToImage( static_cast<DComplexImage&>( *image ) ); break;
2593  }
2594  else
2595  switch ( image.BitsPerSample() )
2596  {
2597  case 8: ToImage( static_cast<UInt8Image&>( *image ) ); break;
2598  case 16: ToImage( static_cast<UInt16Image&>( *image ) ); break;
2599  case 32: ToImage( static_cast<UInt32Image&>( *image ) ); break;
2600  }
2601  }
2602 
2603 #endif // !__PCL_NO_MATRIX_IMAGE_RENDERING
2604 
2605 #ifndef __PCL_NO_MATRIX_IMAGE_CONVERSION
2606 
2644  template <class P>
2645  static GenericMatrix FromImage( const GenericImage<P>& image, const Rect& rect = Rect( 0 ), int channel = -1 )
2646  {
2647  Rect r = rect;
2648  if ( !image.ParseSelection( r, channel ) )
2649  return GenericMatrix();
2650 
2651  if ( r == image.Bounds() )
2652  {
2653  GenericMatrix M( image.Height(), image.Width() );
2654  const typename P::sample* v = image[channel];
2655  block_iterator a = M.m_data->Begin();
2656  const_block_iterator b = M.m_data->End();
2657  for ( ; a < b; ++a, ++v )
2658  P::FromSample( *a, *v );
2659  return M;
2660  }
2661  else
2662  {
2663  int w = r.Width();
2664  int h = r.Height();
2665  GenericMatrix M( h, w );
2666  block_iterator a = M.m_data->Begin();
2667  for ( int i = r.y0; i < r.y1; ++i )
2668  {
2669  const typename P::sample* v = image.PixelAddress( r.x0, i, channel );
2670  for ( int j = 0; j < w; ++j )
2671  P::FromSample( *a++, *v++ );
2672  }
2673  return M;
2674  }
2675  }
2676 
2691  static GenericMatrix FromImage( const ImageVariant& image, const Rect& rect = Rect( 0 ), int channel = -1 )
2692  {
2693  if ( image )
2694  if ( image.IsFloatSample() )
2695  switch ( image.BitsPerSample() )
2696  {
2697  case 32: return FromImage( static_cast<const Image&>( *image ), rect, channel );
2698  case 64: return FromImage( static_cast<const DImage&>( *image ), rect, channel );
2699  }
2700  else if ( image.IsComplexSample() )
2701  switch ( image.BitsPerSample() )
2702  {
2703  case 32: return FromImage( static_cast<const ComplexImage&>( *image ), rect, channel );
2704  case 64: return FromImage( static_cast<const DComplexImage&>( *image ), rect, channel );
2705  }
2706  else
2707  switch ( image.BitsPerSample() )
2708  {
2709  case 8: return FromImage( static_cast<const UInt8Image&>( *image ), rect, channel );
2710  case 16: return FromImage( static_cast<const UInt16Image&>( *image ), rect, channel );
2711  case 32: return FromImage( static_cast<const UInt32Image&>( *image ), rect, channel );
2712  }
2713 
2714  return GenericMatrix();
2715  }
2716 
2717 #endif // !__PCL_NO_MATRIX_IMAGE_CONVERSION
2718 
2719 private:
2720 
2726  struct Data : public ReferenceCounter
2727  {
2728  int n = 0;
2729  int m = 0;
2730  block_iterator* v = nullptr;
2731 
2732  Data() = default;
2733 
2734  Data( int rows, int cols )
2735  {
2736  if ( rows > 0 && cols > 0 )
2737  Allocate( rows, cols );
2738  }
2739 
2740  ~Data()
2741  {
2742  Deallocate();
2743  }
2744 
2745  int Rows() const
2746  {
2747  return n;
2748  }
2749 
2750  int Cols() const
2751  {
2752  return m;
2753  }
2754 
2755  size_type NumberOfElements() const
2756  {
2757  return size_type( n )*size_type( m );
2758  }
2759 
2760  size_type Size() const
2761  {
2762  return NumberOfElements()*sizeof( element );
2763  }
2764 
2765  block_iterator Begin() const
2766  {
2767  return (v != nullptr) ? *v : nullptr;
2768  }
2769 
2770  block_iterator End() const
2771  {
2772  return (v != nullptr) ? *v + NumberOfElements() : nullptr;
2773  }
2774 
2775  element& Element( int i, int j ) const
2776  {
2777  return v[i][j];
2778  }
2779 
2780  void Allocate( int rows, int cols )
2781  {
2782  n = rows;
2783  m = cols;
2784  v = new block_iterator[ n ];
2785  *v = new element[ NumberOfElements() ];
2786  for ( int i = 1; i < n; ++i )
2787  v[i] = v[i-1] + m;
2788  }
2789 
2790  void Deallocate()
2791  {
2792  PCL_PRECONDITION( refCount == 0 )
2793  if ( v != nullptr )
2794  {
2795  delete [] *v;
2796  delete [] v;
2797  v = nullptr;
2798  n = m = 0;
2799  }
2800  }
2801  };
2802 
2807  Data* m_data = nullptr;
2808 
2813  void DetachFromData()
2814  {
2815  if ( !m_data->Detach() )
2816  delete m_data;
2817  }
2818 
2823  static bool Invert( block_iterator* Ai, const_block_iterator const* A, int n )
2824  {
2825  switch ( n )
2826  {
2827  case 1:
2828  if ( 1 + **A == 1 )
2829  return false;
2830  **Ai = 1/(**A);
2831  break;
2832  case 2:
2833  {
2834  const_block_iterator A0 = A[0];
2835  const_block_iterator A1 = A[1];
2836  element d = A0[0]*A1[1] - A0[1]*A1[0];
2837  if ( 1 + d == 1 )
2838  return false;
2839  Ai[0][0] = A1[1]/d;
2840  Ai[0][1] = -A0[1]/d;
2841  Ai[1][0] = -A1[0]/d;
2842  Ai[1][1] = A0[0]/d;
2843  }
2844  break;
2845  case 3:
2846  {
2847  const_block_iterator A0 = A[0];
2848  const_block_iterator A1 = A[1];
2849  const_block_iterator A2 = A[2];
2850  element d1 = A1[1]*A2[2] - A1[2]*A2[1];
2851  element d2 = A1[2]*A2[0] - A1[0]*A2[2];
2852  element d3 = A1[0]*A2[1] - A1[1]*A2[0];
2853  element d = A0[0]*d1 + A0[1]*d2 + A0[2]*d3;
2854  if ( 1 + d == 1 )
2855  return false;
2856  Ai[0][0] = d1/d;
2857  Ai[0][1] = (A2[1]*A0[2] - A2[2]*A0[1])/d;
2858  Ai[0][2] = (A0[1]*A1[2] - A0[2]*A1[1])/d;
2859  Ai[1][0] = d2/d;
2860  Ai[1][1] = (A2[2]*A0[0] - A2[0]*A0[2])/d;
2861  Ai[1][2] = (A0[2]*A1[0] - A0[0]*A1[2])/d;
2862  Ai[2][0] = d3/d;
2863  Ai[2][1] = (A2[0]*A0[1] - A2[1]*A0[0])/d;
2864  Ai[2][2] = (A0[0]*A1[1] - A0[1]*A1[0])/d;
2865  }
2866  break;
2867  default: // ?!
2868  return false;
2869  }
2870  return true;
2871  }
2872 };
2873 
2874 // ### N.B.: Visual C++ is unable to compile the next two member functions if
2875 // the following external function declarations are placed within the
2876 // member function bodies. This forces us to implement them after the
2877 // GenericMatrix<> class declaration.
2878 
2881 
2882 template <typename T> inline
2884 {
2885  if ( Rows() != Cols() || Rows() == 0 )
2886  throw Error( "Invalid matrix inversion: Non-square or empty matrix." );
2887 
2888  /*
2889  * - Use direct formulae to invert 1x1, 2x2 and 3x3 matrices.
2890  * - Use Gauss-Jordan elimination to invert 4x4 and larger matrices.
2891  */
2892  if ( Rows() < 4 )
2893  {
2894  GenericMatrix Ai( Rows(), Rows() );
2895  if ( !Invert( Ai.m_data->v, m_data->v, Rows() ) )
2896  throw Error( "Invalid matrix inversion: Singular matrix." );
2897  return Ai;
2898  }
2899  else
2900  {
2901  GenericMatrix Ai( *this );
2902  GenericMatrix B = UnitMatrix( Rows() );
2903  InPlaceGaussJordan( Ai, B );
2904  return Ai;
2905  }
2906 }
2907 
2908 template <typename T> inline
2910 {
2911  if ( Rows() <= 3 )
2912  Transfer( Inverse() );
2913  else
2914  {
2915  if ( Rows() != Cols() || Rows() == 0 )
2916  throw Error( "Invalid matrix inversion: Non-square or empty matrix." );
2917  EnsureUnique();
2918  GenericMatrix B = UnitMatrix( Rows() );
2919  InPlaceGaussJordan( *this, B );
2920  }
2921 }
2922 
2923 // ----------------------------------------------------------------------------
2924 
2940 template <typename T> inline
2942 {
2943  if ( B.Rows() != A.Rows() || B.Cols() != A.Cols() )
2944  throw Error( "Invalid matrix addition." );
2945  GenericMatrix<T> R( A.Rows(), A.Cols() );
2946  typename GenericMatrix<T>::block_iterator r = R.Begin();
2947  typename GenericMatrix<T>::const_block_iterator s = R.End();
2950  for ( ; r < s; ++r, ++a, ++b )
2951  *r = *a + *b;
2952  return R;
2953 }
2954 
2960 template <typename T> inline
2962 {
2963  GenericMatrix<T> R( A.Rows(), A.Cols() );
2964  typename GenericMatrix<T>::block_iterator r = R.Begin();
2965  typename GenericMatrix<T>::const_block_iterator s = R.End();
2967  for ( ; r < s; ++r, ++a )
2968  *r = *a + x;
2969  return R;
2970 }
2971 
2980 template <typename T> inline
2982 {
2983  return A + x;
2984 }
2985 
2986 // ----------------------------------------------------------------------------
2987 
2996 template <typename T> inline
2998 {
2999  if ( B.Rows() != A.Rows() || B.Cols() != A.Cols() )
3000  throw Error( "Invalid matrix subtraction." );
3001  GenericMatrix<T> R( A.Rows(), A.Cols() );
3002  typename GenericMatrix<T>::block_iterator r = R.Begin();
3003  typename GenericMatrix<T>::const_block_iterator s = R.End();
3006  for ( ; r < s; ++r, ++a, ++b )
3007  *r = *a - *b;
3008  return R;
3009 }
3010 
3016 template <typename T> inline
3018 {
3019  GenericMatrix<T> R( A.Rows(), A.Cols() );
3020  typename GenericMatrix<T>::block_iterator r = R.Begin();
3021  typename GenericMatrix<T>::const_block_iterator s = R.End();
3023  for ( ; r < s; ++r, ++a )
3024  *r = *a - x;
3025  return R;
3026 }
3027 
3037 template <typename T> inline
3039 {
3040  GenericMatrix<T> R( A.Rows(), A.Cols() );
3041  typename GenericMatrix<T>::block_iterator r = R.Begin();
3042  typename GenericMatrix<T>::const_block_iterator s = R.End();
3044  for ( ; r < s; ++r, ++a )
3045  *r = x - *a;
3046  return R;
3047 }
3048 
3049 // ----------------------------------------------------------------------------
3050 
3063 template <typename T> inline
3065 {
3066  int n = A.Rows();
3067  int m = B.Cols();
3068  int p = A.Cols();
3069  if ( B.Rows() != p )
3070  throw Error( "Invalid matrix multiplication." );
3071  GenericMatrix<T> R( n, m );
3072  for ( int i = 0; i < n; ++i )
3073  for ( int j = 0; j < m; ++j )
3074  {
3075  T rij = 0;
3076  for ( int k = 0; k < p; ++k )
3077  rij += A[i][k] * B[k][j];
3078  R[i][j] = rij;
3079  }
3080  return R;
3081 }
3082 
3095 template <typename T> inline
3097 {
3098  int n = A.Cols();
3099  int m = A.Rows();
3100  if ( x.Length() != n )
3101  throw Error( "Invalid matrix-vector multiplication." );
3102  GenericVector<T> r( m );
3103  for ( int i = 0; i < m; ++i )
3104  {
3105  T ri = 0;
3106  for ( int j = 0; j < n; ++j )
3107  ri += A[i][j] * x[j];
3108  r[i] = ri;
3109  }
3110  return r;
3111 }
3112 
3118 template <typename T> inline
3120 {
3121  GenericMatrix<T> R( A.Rows(), A.Cols() );
3122  typename GenericMatrix<T>::block_iterator r = R.Begin();
3123  typename GenericMatrix<T>::const_block_iterator s = R.End();
3125  for ( ; r < s; ++r, ++a )
3126  *r = *a * x;
3127  return R;
3128 }
3129 
3138 template <typename T> inline
3140 {
3141  return A * x;
3142 }
3143 
3144 // ----------------------------------------------------------------------------
3145 
3151 template <typename T> inline
3153 {
3154  GenericMatrix<T> R( A.Rows(), A.Cols() );
3155  typename GenericMatrix<T>::block_iterator r = R.Begin();
3156  typename GenericMatrix<T>::const_block_iterator s = R.End();
3158  for ( ; r < s; ++r, ++a )
3159  *r = *a / x;
3160  return R;
3161 }
3162 
3171 template <typename T> inline
3173 {
3174  GenericMatrix<T> R( A.Rows(), A.Cols() );
3175  typename GenericMatrix<T>::block_iterator r = R.Begin();
3176  typename GenericMatrix<T>::const_block_iterator s = R.End();
3178  for ( ; r < s; ++r, ++a )
3179  *r = x / *a;
3180  return R;
3181 }
3182 
3183 // ----------------------------------------------------------------------------
3184 
3190 template <typename T> inline
3192 {
3193  GenericMatrix<T> R( A.Rows(), A.Cols() );
3194  typename GenericMatrix<T>::block_iterator r = R.Begin();
3195  typename GenericMatrix<T>::const_block_iterator s = R.End();
3197  for ( ; r < s; ++r, ++a )
3198  *r = pcl::Pow( *a, x );
3199  return R;
3200 }
3201 
3210 template <typename T> inline
3212 {
3213  GenericMatrix<T> R( A.Rows(), A.Cols() );
3214  typename GenericMatrix<T>::block_iterator r = R.Begin();
3215  typename GenericMatrix<T>::const_block_iterator s = R.End();
3217  for ( ; r < s; ++r, ++a )
3218  *r = pcl::Pow( x, *a );
3219  return R;
3220 }
3221 
3222 // ----------------------------------------------------------------------------
3223 
3224 #ifndef __PCL_NO_MATRIX_INSTANTIATE
3225 
3238 
3247 typedef I8Matrix CharMatrix;
3248 
3257 
3266 typedef UI8Matrix ByteMatrix;
3267 
3276 
3285 
3294 
3303 typedef I32Matrix IMatrix;
3304 
3313 
3322 typedef UI32Matrix UIMatrix;
3323 
3332 
3341 
3350 
3359 typedef F32Matrix FMatrix;
3360 
3369 
3378 typedef F64Matrix DMatrix;
3379 
3388 typedef DMatrix Matrix;
3389 
3398 
3407 
3408 #ifndef _MSC_VER
3409 
3422 
3434 typedef F80Matrix LDMatrix;
3435 
3436 #endif // !_MSC_VER
3437 
3438 #endif // !__PCL_NO_MATRIX_INSTANTIATE
3439 
3440 // ----------------------------------------------------------------------------
3441 
3442 } // pcl
3443 
3444 #endif // __PCL_Matrix_h
3445 
3446 // ----------------------------------------------------------------------------
3447 // EOF pcl/Matrix.h - Released 2019-11-07T10:59:34Z
void ReverseSort()
Definition: Matrix.h:1799
static GenericMatrix CrossPowerSpectrumMatrix(const GenericMatrix &A, const GenericMatrix &B)
Definition: Matrix.h:2520
element MaxElement(Point &p) const
Definition: Matrix.h:1928
vector ColVector(int j) const
Definition: Matrix.h:1297
double Mean(const T *i, const T *j)
Definition: Math.h:2351
void CrossPowerSpectrumMatrix(Complex< T > *i, const Complex< T > *j, const Complex< T > *a, const Complex< T > *b)
Definition: Complex.h:1005
uint64 Hash64(const void *data, size_type size, uint64 seed=0)
Definition: Math.h:3551
size_type Size() const
Definition: Matrix.h:984
uint64 Hash(uint64 seed=0) const
Definition: Matrix.h:2377
GenericMatrix(GenericMatrix &&x)
Definition: Matrix.h:246
GenericVector< element > vector
Definition: Matrix.h:134
void RotateY(double sphi, double cphi)
Definition: Matrix.h:1552
Complex< T1 > operator*(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:539
GenericImage< P > operator^(const GenericImage< P > &image, T scalar)
Definition: Image.h:16199
void Transfer(GenericMatrix &x)
Definition: Matrix.h:440
double BiweightMidvariance(double center, double sigma, int k=9) const
Definition: Matrix.h:2211
BI LinearSearchLast(BI i, BI j, const T &v)
Definition: Search.h:129
Complex< T1 > operator+(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:455
void Transfer(GenericMatrix &&x)
Definition: Matrix.h:459
bool IsFloatSample() const
Definition: ImageVariant.h:435
GenericMatrix Sorted() const
Definition: Matrix.h:1789
double StableAvgDev(const T *i, const T *j, double center)
Definition: Math.h:2770
void Swap(GenericPoint< T > &p1, GenericPoint< T > &p2)
Definition: Point.h:1402
virtual ~GenericMatrix()
Definition: Matrix.h:357
bool ParseSelection(Rect &rect, int &firstChannel, int &lastChannel) const
uint32 Hash32(uint32 seed=0) const
Definition: Matrix.h:2368
double Mean() const
Definition: Matrix.h:2010
double Median() const
Definition: Matrix.h:2077
GenericMatrix Sqr() const
Definition: Matrix.h:741
GenericMatrix(const GenericImage< P > &image, const Rect &rect=Rect(0), int channel=-1)
Definition: Matrix.h:324
vector ColumnVector(int j) const
Definition: Matrix.h:1284
static GenericMatrix FromImage(const ImageVariant &image, const Rect &rect=Rect(0), int channel=-1)
Definition: Matrix.h:2691
void Truncate(const element &f0=element(0), const element &f1=element(1))
Definition: Matrix.h:1701
double Median(T *i, T *j)
Definition: Math.h:2514
vector RowVector(int i) const
Definition: Matrix.h:1273
int Compare(FI1 i1, FI1 j1, FI2 i2, FI2 j2)
Definition: Utility.h:639
int BitsPerSample() const
Definition: ImageVariant.h:452
PCL root namespace.
Definition: AbstractImage.h:76
void SetSqrt()
Definition: Matrix.h:802
GenericMatrix(const GenericMatrix &x)
Definition: Matrix.h:238
GenericMatrix Transpose() const
Definition: Matrix.h:1392
double Qn(T *x, T *xn)
Definition: Math.h:3186
32-bit signed integer matrix.
Rect Bounds() const
bool SameDimensions(const GenericMatrix &x) const
Definition: Matrix.h:1019
32-bit floating point real matrix.
64-bit integer matrix.
const_block_iterator Begin() const
Definition: Matrix.h:1122
32-bit floating point complex matrix.
64-bit floating point real matrix.
64-bit floating point complex matrix.
bool SameElements(const GenericMatrix &x) const
Definition: Matrix.h:1033
GenericMatrix(int rows, int cols)
Definition: Matrix.h:168
FI LinearSearch(FI i, FI j, const T &v)
Definition: Search.h:91
double Sum() const
Definition: Matrix.h:1947
int Rows() const
Definition: Matrix.h:907
const_block_iterator end() const
Definition: Matrix.h:1264
int Length() const
Definition: Vector.h:1502
const_block_iterator End() const
Definition: Matrix.h:1185
8-bit signed integer matrix.
Complex< T > Sqrt(const Complex< T > &c)
Definition: Complex.h:665
FI MaxItem(FI i, FI j)
Definition: Utility.h:479
element MinElement() const
Definition: Matrix.h:1884
double BiweightMidvariance(const T *x, const T *xn, double center, double sigma, int k=9)
Definition: Math.h:3343
64-bit floating point real matrix.
void Sort(BI i, BI j)
Definition: Sort.h:534
bool IsValid() const
Definition: Matrix.h:947
double BendMidvariance(double center, double beta=0.2) const
Definition: Matrix.h:2263
Implements a generic, two-dimensional, shared or local image.
Definition: Image.h:275
GenericMatrix RotatedY(double phi) const
Definition: Matrix.h:1602
int Width() const
Definition: ImageGeometry.h:90
size_type NumberOfElements() const
Definition: Matrix.h:975
double Sum(const T *i, const T *j)
Definition: Math.h:2215
bool Equal(FI1 i1, FI2 i2, FI2 j2)
Definition: Utility.h:592
int Height() const
Definition: ImageGeometry.h:98
void RotateY(double phi)
Definition: Matrix.h:1579
double StableSum() const
Definition: Matrix.h:1958
GenericImage & AllocateData(int width, int height, int numberOfChannels=1, color_space colorSpace=ColorSpace::Gray)
Definition: Image.h:6535
block_iterator end()
Definition: Matrix.h:1256
Unsigned integer matrix.
double Median()
Definition: Matrix.h:2059
GenericMatrix Flipped() const
Definition: Matrix.h:1442
void SetCol(int j, const V &c)
Definition: Matrix.h:1333
block_iterator RowPtr(int i)
Definition: Matrix.h:1231
GenericMatrix Sorted(BP p) const
Definition: Matrix.h:1833
Generic vector of arbitrary length.
Definition: Vector.h:106
double StableMean(const T *i, const T *j)
Definition: Math.h:2370
static GenericMatrix FromImage(const GenericImage< P > &image, const Rect &rect=Rect(0), int channel=-1)
Definition: Matrix.h:2645
constexpr const T & Range(const T &x, const T &a, const T &b)
Definition: Utility.h:190
GenericMatrix Inverse() const
Definition: Matrix.h:2883
size_t size_type
Definition: Defs.h:543
double Variance(const T *i, const T *j, double center)
Definition: Math.h:2395
bool IsUnique() const
Definition: Matrix.h:871
8-bit unsigned integer matrix.
const_block_iterator FindFirst(const element &x) const
Definition: Matrix.h:1855
Acts like a union for all types of images in PCL, with optional class-wide ownership of transported i...
Definition: ImageVariant.h:317
GenericMatrix(const T1 &a00, const T1 &a01, const T1 &a02, const T1 &a10, const T1 &a11, const T1 &a12, const T1 &a20, const T1 &a21, const T1 &a22)
Definition: Matrix.h:223
64-bit floating point real matrix.
T Abs(const Complex< T > &c)
Definition: Complex.h:420
double Modulus() const
Definition: Matrix.h:1968
void Assign(const GenericMatrix &x)
Definition: Matrix.h:412
double StableSum(const T *i, const T *j)
Definition: Math.h:2234
int Cols() const
Definition: Matrix.h:916
unsigned long long uint64
Definition: Defs.h:616
static GenericMatrix FromColumnVector(const vector &c)
Definition: Matrix.h:1356
double MAD(const T *i, const T *j, double center)
Definition: Math.h:2865
sample * PixelAddress(int x, int y, int channel=0)
Definition: Image.h:6915
GenericMatrix Sqrt() const
Definition: Matrix.h:785
static GenericMatrix PhaseCorrelationMatrix(const GenericMatrix &A, const GenericMatrix &B)
Definition: Matrix.h:2507
64-bit unsigned integer matrix.
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
const_block_iterator FindLast(const element &x) const
Definition: Matrix.h:1864
void SetColumn(int j, const V &c)
Definition: Matrix.h:1319
static GenericMatrix FromColVector(const vector &c)
Definition: Matrix.h:1367
GenericMatrix RotatedZ(double sphi, double cphi) const
Definition: Matrix.h:1671
8-bit signed integer matrix.
GenericMatrix Rescaled(const element &f0=element(0), const element &f1=element(1)) const
Definition: Matrix.h:1770
16-bit signed integer matrix.
double StableSumOfSquares(const T *i, const T *j)
Definition: Math.h:2325
bool IsComplexSample() const
Definition: ImageVariant.h:443
constexpr const T & Min(const T &a, const T &b)
Definition: Utility.h:90
static GenericMatrix FromRowVector(const vector &r)
Definition: Matrix.h:1343
block_iterator * DataPtr()
Definition: Matrix.h:1216
bool IsAliasOf(const GenericMatrix &x) const
Definition: Matrix.h:880
void PCL_FUNC InPlaceGaussJordan(Matrix &A, Matrix &B)
block_iterator Begin()
Definition: Matrix.h:1107
Complex< T1 > operator-(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:495
element MinElement(Point &p) const
Definition: Matrix.h:1897
void ToImage(ImageVariant &image) const
Definition: Matrix.h:2577
void RotateX(double sphi, double cphi)
Definition: Matrix.h:1472
element & Element(int i, int j)
Definition: Matrix.h:1054
32-bit unsigned integer matrix.
void Rescale(const element &f0=element(0), const element &f1=element(1))
Definition: Matrix.h:1738
uint32 Hash32(const void *data, size_type size, uint32 seed=0)
Definition: Math.h:3699
double StableAvgDev(double center) const
Definition: Matrix.h:2116
element * block_iterator
Definition: Matrix.h:141
bool operator<(const Array< T, A > &x1, const Array< T, A > &x2)
Definition: Array.h:2086
Generic dynamic matrix of arbitrary dimensions.
Definition: Matrix.h:122
ImageVariant & CreateFloatImage(int bitSize=32)
A simple exception with an associated error message.
Definition: Exception.h:213
double StableModulus() const
Definition: Matrix.h:1979
GenericMatrix(const GenericMatrix &x, int i0, int j0, int rows, int cols)
Definition: Matrix.h:273
void ToImage(GenericImage< P > &image) const
Definition: Matrix.h:2552
GenericMatrix RotatedY(double sphi, double cphi) const
Definition: Matrix.h:1591
const_block_iterator ConstEnd() const
Definition: Matrix.h:1193
void SetRow(int i, const V &r)
Definition: Matrix.h:1307
void EnsureUnique()
Definition: Matrix.h:892
double BendMidvariance(const T *x, const T *xn, double center, double beta=0.2)
Definition: Math.h:3402
double StableMean() const
Definition: Matrix.h:2021
double SumOfSquares() const
Definition: Matrix.h:1989
block_iterator begin()
Definition: Matrix.h:1240
Root base class of all PCL containers of objects.
Definition: Container.h:77
double Variance() const
Definition: Matrix.h:2031
GenericMatrix(const T1 *a, int rows, int cols)
Definition: Matrix.h:200
const_block_iterator Find(const element &x) const
Definition: Matrix.h:1844
static GenericMatrix UnitMatrix(int n)
Definition: Matrix.h:1378
32-bit signed integer matrix.
GenericMatrix ReverseSorted() const
Definition: Matrix.h:1809
double MAD() const
Definition: Matrix.h:2182
element MaxElement() const
Definition: Matrix.h:1915
double Sn(T *x, T *xn)
Definition: Math.h:2939
const_block_iterator ConstBegin() const
Definition: Matrix.h:1130
FI MinItem(FI i, FI j)
Definition: Utility.h:441
void SinCos(T x, T &sx, T &cx)
Definition: Math.h:938
16-bit unsigned integer matrix.
void PhaseCorrelationMatrix(Complex< T > *i, const Complex< T > *j, const Complex< T > *a, const Complex< T > *b)
Definition: Complex.h:974
void RotateX(double phi)
Definition: Matrix.h:1499
bool operator==(const Array< T, A > &x1, const Array< T, A > &x2)
Definition: Array.h:2075
Complex< T1 > operator/(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:583
32-bit integer point on the plane.
bool Contains(const element &x) const
Definition: Matrix.h:1873
32-bit floating point real matrix.
double Modulus(const T *i, const T *j)
Definition: Math.h:2259
double BiweightMidvariance(int k=9) const
Definition: Matrix.h:2234
double Sn() const
Definition: Matrix.h:2313
int Columns() const
Definition: Matrix.h:927
block_iterator End()
Definition: Matrix.h:1169
GenericMatrix Abs() const
Definition: Matrix.h:829
ptrdiff_t distance_type
Definition: Defs.h:549
uint64 Hash64(uint64 seed=0) const
Definition: Matrix.h:2355
friend void Swap(GenericMatrix &x1, GenericMatrix &x2)
Definition: Matrix.h:472
void RotateZ(double phi)
Definition: Matrix.h:1659
double StdDev() const
Definition: Matrix.h:2042
80-bit extended precision floating point real matrix.
GenericMatrix RotatedZ(double phi) const
Definition: Matrix.h:1682
double StableSumOfSquares() const
Definition: Matrix.h:2000
double AvgDev() const
Definition: Matrix.h:2133
void RotateZ(double sphi, double cphi)
Definition: Matrix.h:1632
double SumOfSquares(const T *i, const T *j)
Definition: Math.h:2303
const_block_iterator begin() const
Definition: Matrix.h:1248
const element & Element(int i, int j) const
Definition: Matrix.h:1064
GenericMatrix RotatedX(double phi) const
Definition: Matrix.h:1522
GenericMatrix Truncated(const element &f0=element(0), const element &f1=element(1)) const
Definition: Matrix.h:1716
double MAD(double center) const
Definition: Matrix.h:2167
double BendMidvariance(double beta=0.2) const
Definition: Matrix.h:2287
double StableAvgDev() const
Definition: Matrix.h:2151
GenericMatrix(const ImageVariant &image, const Rect &rect=Rect(0), int channel=-1)
Definition: Matrix.h:344
double StableModulus(const T *i, const T *j)
Definition: Math.h:2278
80-bit extended precision floating point real matrix.
unsigned int uint32
Definition: Defs.h:600
void Sort(BP p)
Definition: Matrix.h:1822
32-bit integer rectangle on the plane.
bool IsEmpty() const
Definition: Matrix.h:956
GenericMatrix RotatedX(double sphi, double cphi) const
Definition: Matrix.h:1511
Thread-safe reference counter for copy-on-write data structures.
Complex< T1 > Pow(const Complex< T1 > &c, T2 x)
Definition: Complex.h:738
double Qn() const
Definition: Matrix.h:2339
GenericMatrix(const element &x, int rows, int cols)
Definition: Matrix.h:180
const element * const_block_iterator
Definition: Matrix.h:148
double StdDev(const T *i, const T *j, double center)
Definition: Math.h:2453
double AvgDev(const T *i, const T *j, double center)
Definition: Math.h:2739
8-bit unsigned integer matrix.
double AvgDev(double center) const
Definition: Matrix.h:2096