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_hxx
19 #define itkVariableSizeMatrix_hxx
20
21 #include "itkVariableSizeMatrix.h"
22 #include "itkNumericTraits.h"
23
24 namespace itk
25 {
26 template< typename T >
27 VariableSizeMatrix< T >
VariableSizeMatrix(unsigned int rows,unsigned int cols)28 ::VariableSizeMatrix(unsigned int rows, unsigned int cols):
29 m_Matrix(rows, cols) {}
30
31 /**
32 * Product by a Vector
33 */
34 template< typename T >
35 Array< T >
36 VariableSizeMatrix< T >
operator *(const Array<T> & vect) const37 ::operator*(const Array< T > & vect) const
38 {
39 unsigned int rows = this->Rows();
40 unsigned int cols = this->Cols();
41
42 if ( vect.Size() != cols )
43 {
44 itkGenericExceptionMacro( << "Matrix with " << this->Cols() << " columns cannot be "
45 << "multiplied with array of length: " << vect.Size() );
46 }
47
48 Array< T > result(rows);
49 for ( unsigned int r = 0; r < rows; r++ )
50 {
51 T sum = NumericTraits< T >::ZeroValue();
52 for ( unsigned int c = 0; c < cols; c++ )
53 {
54 sum += m_Matrix(r, c) * vect[c];
55 }
56 result[r] = sum;
57 }
58 return result;
59 }
60
61 /**
62 * Product by a matrix
63 */
64 template< typename T >
65 VariableSizeMatrix< T >
66 VariableSizeMatrix< T >
operator *(const Self & matrix) const67 ::operator*(const Self & matrix) const
68 {
69 if ( this->Cols() != matrix.Rows() )
70 {
71 itkGenericExceptionMacro(<< "Matrix with size ("
72 << this->Rows() << ","
73 << this->Cols()
74 << ") cannot be multiplied by matrix with size ("
75 << matrix.Rows() << "," << matrix.Cols() << ")");
76 }
77 Self result;
78 result = m_Matrix * matrix.m_Matrix;
79 return result;
80 }
81
82 /**
83 * Matrix Addition
84 */
85 template< typename T >
86 VariableSizeMatrix< T >
87 VariableSizeMatrix< T >
operator +(const Self & matrix) const88 ::operator+(const Self & matrix) const
89 {
90 if ( ( matrix.Rows() != this->Rows() )
91 || ( matrix.Cols() != this->Cols() ) )
92 {
93 itkGenericExceptionMacro(<< "Matrix with size (" << matrix.Rows() << ","
94 << matrix.Cols() << ") cannot be added to a matrix with size ("
95 << this->Rows() << "," << this->Cols() << ")");
96 }
97
98 Self result( this->Rows(), this->Cols() );
99 for ( unsigned int r = 0; r < this->Rows(); r++ )
100 {
101 for ( unsigned int c = 0; c < this->Cols(); c++ )
102 {
103 result.m_Matrix(r, c) = m_Matrix(r, c) + matrix.m_Matrix(r, c);
104 }
105 }
106 return result;
107 }
108
109 /**
110 * Matrix Addition in-place
111 */
112 template< typename T >
113 const VariableSizeMatrix< T > &
114 VariableSizeMatrix< T >
operator +=(const Self & matrix)115 ::operator+=(const Self & matrix)
116 {
117 if ( ( matrix.Rows() != this->Rows() )
118 || ( matrix.Cols() != this->Cols() ) )
119 {
120 itkGenericExceptionMacro(<< "Matrix with size (" << matrix.Rows() << ","
121 << matrix.Cols() << ") cannot be added to a matrix with size ("
122 << this->Rows() << "," << this->Cols() << ")");
123 }
124
125 for ( unsigned int r = 0; r < this->Rows(); r++ )
126 {
127 for ( unsigned int c = 0; c < this->Cols(); c++ )
128 {
129 m_Matrix(r, c) += matrix.m_Matrix(r, c);
130 }
131 }
132 return *this;
133 }
134
135 /**
136 * Matrix Subtraction
137 */
138 template< typename T >
139 VariableSizeMatrix< T >
140 VariableSizeMatrix< T >
operator -(const Self & matrix) const141 ::operator-(const Self & matrix) const
142 {
143 if ( ( matrix.Rows() != this->Rows() )
144 || ( matrix.Cols() != this->Cols() ) )
145 {
146 itkGenericExceptionMacro(<< "Matrix with size (" << matrix.Rows() << ","
147 << matrix.Cols() << ") cannot be subtracted from matrix with size ("
148 << this->Rows() << "," << this->Cols() << ")");
149 }
150
151 Self result( this->Rows(), this->Cols() );
152 for ( unsigned int r = 0; r < this->Rows(); r++ )
153 {
154 for ( unsigned int c = 0; c < this->Cols(); c++ )
155 {
156 result.m_Matrix(r, c) = m_Matrix(r, c) - matrix.m_Matrix(r, c);
157 }
158 }
159 return result;
160 }
161
162 /**
163 * Matrix subtraction in-place
164 */
165 template< typename T >
166 const VariableSizeMatrix< T > &
167 VariableSizeMatrix< T >
operator -=(const Self & matrix)168 ::operator-=(const Self & matrix)
169 {
170 if ( ( matrix.Rows() != this->Rows() )
171 || ( matrix.Cols() != this->Cols() ) )
172 {
173 itkGenericExceptionMacro(<< "Matrix with size (" << matrix.Rows() << ","
174 << matrix.Cols() << ") cannot be subtracted from matrix with size ("
175 << this->Rows() << "," << this->Cols() << ")");
176 }
177
178 for ( unsigned int r = 0; r < this->Rows(); r++ )
179 {
180 for ( unsigned int c = 0; c < this->Cols(); c++ )
181 {
182 m_Matrix(r, c) -= matrix.m_Matrix(r, c);
183 }
184 }
185 return *this;
186 }
187
188 template< typename T >
189 VariableSizeMatrix< T > &
190 VariableSizeMatrix< T >
operator -()191 ::operator-()
192 {
193 for ( unsigned int r = 0; r < this->Rows(); r++ )
194 {
195 for ( unsigned int c = 0; c < this->Cols(); c++ )
196 {
197 m_Matrix(r, c) = -m_Matrix(r, c);
198 }
199 }
200 return *this;
201 }
202
203 /**
204 * Product by a vnl_matrix
205 */
206 template< typename T >
207 vnl_matrix< T >
208 VariableSizeMatrix< T >
operator *(const vnl_matrix<T> & matrix) const209 ::operator*(const vnl_matrix< T > & matrix) const
210 {
211 return m_Matrix * matrix;
212 }
213
214 /**
215 * Product by a matrix
216 */
217 template< typename T >
218 void
219 VariableSizeMatrix< T >
operator *=(const Self & matrix)220 ::operator*=(const Self & matrix)
221 {
222 m_Matrix *= matrix.m_Matrix;
223 }
224
225 /**
226 * Product by a vnl_matrix
227 */
228 template< typename T >
229 void
230 VariableSizeMatrix< T >
operator *=(const vnl_matrix<T> & matrix)231 ::operator*=(const vnl_matrix< T > & matrix)
232 {
233 m_Matrix *= matrix;
234 }
235
236 /**
237 * Product by a vnl_vector
238 */
239 template< typename T >
240 vnl_vector< T >
241 VariableSizeMatrix< T >
operator *(const vnl_vector<T> & vc) const242 ::operator*(const vnl_vector< T > & vc) const
243 {
244 return m_Matrix * vc;
245 }
246 } // end namespace itk
247
248 #endif
249