1 /*=========================================================================
2  *
3  *  Copyright Insight Software Consortium
4  *
5  *  Licensed under the Apache License, Version 2.0 (the "License");
6  *  you may not use this file except in compliance with the License.
7  *  You may obtain a copy of the License at
8  *
9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  *  Unless required by applicable law or agreed to in writing, software
12  *  distributed under the License is distributed on an "AS IS" BASIS,
13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  *  See the License for the specific language governing permissions and
15  *  limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkVariableSizeMatrix_h
19 #define itkVariableSizeMatrix_h
20 
21 #include "itkPoint.h"
22 #include "itkCovariantVector.h"
23 #include "vnl/vnl_matrix_fixed.h"
24 #include "vnl/algo/vnl_matrix_inverse.h"
25 #include "vnl/vnl_transpose.h"
26 #include "vnl/vnl_matrix.h"
27 #include "itkArray.h"
28 #include "itkMath.h"
29 
30 namespace itk
31 {
32 /** \class VariableSizeMatrix
33  * \brief A templated class holding a M x N size Matrix.
34  *
35  * This class contains a vnl_matrix in order
36  * to make all the vnl mathematical methods available. This class is meant to be
37  * used when the matrix length cannot be determined at compile time.
38  *
39  * \ingroup DataRepresentation
40  * \ingroup ITKCommon
41  */
42 
43 template< typename T >
44 class ITK_TEMPLATE_EXPORT VariableSizeMatrix
45 {
46 public:
47   /** Standard class type aliases. */
48   using Self = VariableSizeMatrix;
49 
50   /** Component value type */
51   using ValueType = T;
52   using ComponentType = T;
53 
54   /** Internal matrix type */
55   using InternalMatrixType = vnl_matrix< T >;
56 
57   /** Matrix by Vector multiplication.  */
58   Array< T > operator *(const Array< T > & vector) const;
59 
60   /** Matrix by Matrix multiplication.  */
61   Self operator *(const Self & matrix) const;
62 
63   /** Matrix addition.  */
64   Self operator+(const Self & matrix) const;
65 
66   const Self & operator+=(const Self & matrix);
67 
68   /** Matrix addition.  */
69   Self operator-(const Self & matrix) const;
70 
71   const Self & operator-=(const Self & matrix);
72 
73   /** negation operator */
74   Self & operator-();
75 
76   /** Matrix by vnl_matrix multiplication.  */
77   vnl_matrix< T > operator *(const vnl_matrix< T > & matrix) const;
78 
79   /** Matrix by Matrix multiplication.  */
80   void operator*=(const Self & matrix);
81 
82   /** Matrix by vnl_matrix multiplication.  */
83   void operator*=(const vnl_matrix< T > & matrix);
84 
85   /** Matrix by vnl_vector multiplication.  */
86   vnl_vector< T > operator *(const vnl_vector< T > & matrix) const;
87 
88   /** Matrix by scalar multiplication.  */
89   void operator*=(const T & value) { m_Matrix *= value; }
90 
91   /** Matrix by scalar multiplication.  */
92   Self operator*(const T & value)
93   {
94     Self result(*this);
95 
96     result *= value;
97     return result;
98   }
99 
100   /** Matrix by scalar division.  */
101   void operator/=(const T & value) { m_Matrix /= value; }
102 
103   /** Matrix by scalar division.  */
104   Self operator/(const T & value)
105   {
106     Self result(*this);
107 
108     result /= value;
109     return result;
110   }
111 
112   /** Return an element of the matrix. */
operator()113   inline T & operator()(unsigned int row, unsigned int col)
114   {
115     return m_Matrix(row, col);
116   }
117 
118   /** Return an element of the matrix. */
operator()119   inline const T & operator()(unsigned int row, unsigned int col) const
120   {
121     return m_Matrix(row, col);
122   }
123 
124   /** Return a row of the matrix. */
125   inline T * operator[](unsigned int i)
126   {
127     return m_Matrix[i];
128   }
129 
130   /** Return a row of the matrix. */
131   inline const T * operator[](unsigned int i) const
132   {
133     return m_Matrix[i];
134   }
135 
136   /** Return the matrix. */
GetVnlMatrix()137   inline InternalMatrixType & GetVnlMatrix()
138   {
139     return m_Matrix;
140   }
141 
142   /** Return the matrix. */
GetVnlMatrix()143   inline const InternalMatrixType & GetVnlMatrix() const
144   {
145     return m_Matrix;
146   }
147 
148   /** Set the matrix to identity. */
SetIdentity()149   inline void SetIdentity()
150   {
151     m_Matrix.set_identity();
152   }
153 
154   /** Fill the matrix with a value. */
Fill(const T & value)155   inline void Fill(const T & value)
156   {
157     m_Matrix.fill(value);
158   }
159 
160   /** Assignment operator. */
161   inline const Self & operator=(const vnl_matrix< T > & matrix)
162   {
163     m_Matrix = matrix;
164     return *this;
165   }
166 
167   /** Comparison operators. */
168   inline bool operator==(const Self & matrix) const;
169 
170   inline bool operator!=(const Self & matrix) const
171   {
172     return !this->operator==(matrix);
173   }
174 
175   /** Assignment operator. */
176   inline const Self & operator=(const Self & matrix)
177   {
178     m_Matrix = matrix.m_Matrix;
179     return *this;
180   }
181 
182   /** Return the inverse matrix. */
GetInverse()183   inline vnl_matrix< T > GetInverse() const
184   {
185     vnl_matrix< T > temp = vnl_matrix_inverse< T >(m_Matrix);
186     return temp;
187   }
188 
189   /** Return the transposed matrix. */
GetTranspose()190   inline vnl_matrix< T > GetTranspose() const
191   {
192     return m_Matrix.transpose();
193   }
194 
195   /** Default constructor. */
VariableSizeMatrix()196   VariableSizeMatrix():m_Matrix() {}
197 
198   VariableSizeMatrix(unsigned int rows, unsigned int cols);
199 
200   /** Copy constructor. */
VariableSizeMatrix(const Self & matrix)201   VariableSizeMatrix(const Self & matrix):m_Matrix(matrix.m_Matrix) {}
202 
203   /** Return number of rows in the matrix */
Rows()204   inline unsigned int Rows() const { return m_Matrix.rows(); }
205 
206   /** Return number of columns in the matrix */
Cols()207   inline unsigned int Cols() const { return m_Matrix.cols(); }
208 
209   /** Set the matrix size. Old data lost. Returns true if size changed. */
SetSize(unsigned int r,unsigned int c)210   inline bool SetSize(unsigned int r, unsigned int c) { return m_Matrix.set_size(r, c); }
211 
212 private:
213   InternalMatrixType m_Matrix;
214 };
215 
216 template< typename T >
217 std::ostream & operator<<(std::ostream & os,
218                                      const VariableSizeMatrix< T > & v)
219 {
220   os << v.GetVnlMatrix(); return os;
221 }
222 
223 /**
224  *  Comparison
225  */
226 template< typename T >
227 inline
228 bool
229 VariableSizeMatrix< T >
230 ::operator==(const Self & matrix) const
231 {
232   if ( ( matrix.Rows() != this->Rows() )
233        || ( matrix.Cols() != this->Cols() ) )
234     {
235     return false;
236     }
237   bool equal = true;
238 
239   for ( unsigned int r = 0; r < this->Rows(); r++ )
240     {
241     for ( unsigned int c = 0; c < this->Cols(); c++ )
242       {
243       if ( Math::NotExactlyEquals(m_Matrix(r, c), matrix.m_Matrix(r, c)) )
244         {
245         equal = false;
246         break;
247         }
248       }
249     }
250   return equal;
251 }
252 } // end namespace itk
253 
254 #ifndef ITK_MANUAL_INSTANTIATION
255 #include "itkVariableSizeMatrix.hxx"
256 #endif
257 
258 #endif
259