go home Home | Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Data Structures | File List | Namespace Members | Data Fields | Globals | Related Pages
vnl_adjugate_fixed.h
Go to the documentation of this file.
1 // This is core/vnl/vnl_adjugate_fixed.h
2 #ifndef vnl_adjugate_fixed_h_
3 #define vnl_adjugate_fixed_h_
4 //:
5 // \file
6 // \brief Calculates adjugate or cofactor matrix of a small vnl_matrix_fixed.
7 // Code is only a small adaptation from Peter Vanroose's vnl_inverse.h
8 // adjugate == inverse/det
9 // cofactor == adjugate^T
10 // \author Floris Berendsen
11 // \date 18 April 2013
12 //
13 // \verbatim
14 // Code is only a small adaptation from Peter Vanroose's vnl_inverse.h
15 // \endverbatim
16 
17 #include <vnl/vnl_matrix_fixed.h>
18 #include <vnl/vnl_vector_fixed.h>
19 #include <vnl/vnl_matrix.h>
20 #include <vnl/vnl_det.h>
21 #include <vcl_cassert.h>
22 
23 //: Calculates adjugate of a small vnl_matrix_fixed (not using svd)
24 // This allows you to write e.g.
25 //
26 // x = vnl_adjugate(A) * b;
27 //
28 // Note that this function is inlined (except for the call to vnl_det()),
29 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
30 // since that one is using svd.
31 //
32 // \relatesalso vnl_matrix_fixed
33 
34 template <class T>
35 vnl_matrix_fixed<T,1,1> vnl_adjugate(vnl_matrix_fixed<T,1,1> const& m)
36 {
37  return vnl_matrix_fixed<T,1,1>(m(0,0));
38 }
39 
40 //: Calculates adjugate of a small vnl_matrix_fixed (not using svd)
41 // This allows you to write e.g.
42 //
43 // x = vnl_adjugate(A) * b;
44 //
45 // Note that this function is inlined (except for the call to vnl_det()),
46 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
47 // since that one is using svd.
48 //
49 // \relatesalso vnl_matrix_fixed
50 
51 template <class T>
52 vnl_matrix_fixed<T,2,2> vnl_adjugate(vnl_matrix_fixed<T,2,2> const& m)
53 {
54  T d[4];
55  d[0] = m(1,1); d[1] = - m(0,1);
56  d[3] = m(0,0); d[2] = - m(1,0);
57  return vnl_matrix_fixed<T,2,2>(d);
58 }
59 
60 //: Calculates adjugate of a small vnl_matrix_fixed (not using svd)
61 // This allows you to write e.g.
62 //
63 // x = vnl_adjugate_fixed(A) * b;
64 //
65 // Note that this function is inlined (except for the call to vnl_det()),
66 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
67 // since that one is using svd.
68 //
69 // \relatesalso vnl_matrix_fixed
70 
71 template <class T>
72 vnl_matrix_fixed<T,3,3> vnl_adjugate(vnl_matrix_fixed<T,3,3> const& m)
73 {
74  T d[9];
75  d[0] = (m(1,1)*m(2,2)-m(1,2)*m(2,1));
76  d[1] = (m(2,1)*m(0,2)-m(2,2)*m(0,1));
77  d[2] = (m(0,1)*m(1,2)-m(0,2)*m(1,1));
78  d[3] = (m(1,2)*m(2,0)-m(1,0)*m(2,2));
79  d[4] = (m(0,0)*m(2,2)-m(0,2)*m(2,0));
80  d[5] = (m(1,0)*m(0,2)-m(1,2)*m(0,0));
81  d[6] = (m(1,0)*m(2,1)-m(1,1)*m(2,0));
82  d[7] = (m(0,1)*m(2,0)-m(0,0)*m(2,1));
83  d[8] = (m(0,0)*m(1,1)-m(0,1)*m(1,0));
84  return vnl_matrix_fixed<T,3,3>(d);
85 }
86 
87 //: Calculates adjugate_fixed of a small vnl_matrix_fixed (not using svd)
88 // This allows you to write e.g.
89 //
90 // x = vnl_adjugate_fixed(A) * b;
91 //
92 // Note that this function is inlined (except for the call to vnl_det()),
93 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
94 // since that one is using svd.
95 //
96 // \relatesalso vnl_matrix_fixed
97 
98 template <class T>
99 vnl_matrix_fixed<T,4,4> vnl_adjugate(vnl_matrix_fixed<T,4,4> const& m)
100 {
101  T d[16];
102  d[0] = m(1,1)*m(2,2)*m(3,3) - m(1,1)*m(2,3)*m(3,2) - m(2,1)*m(1,2)*m(3,3)
103  + m(2,1)*m(1,3)*m(3,2) + m(3,1)*m(1,2)*m(2,3) - m(3,1)*m(1,3)*m(2,2);
104  d[1] = -m(0,1)*m(2,2)*m(3,3) + m(0,1)*m(2,3)*m(3,2) + m(2,1)*m(0,2)*m(3,3)
105  - m(2,1)*m(0,3)*m(3,2) - m(3,1)*m(0,2)*m(2,3) + m(3,1)*m(0,3)*m(2,2);
106  d[2] = m(0,1)*m(1,2)*m(3,3) - m(0,1)*m(1,3)*m(3,2) - m(1,1)*m(0,2)*m(3,3)
107  + m(1,1)*m(0,3)*m(3,2) + m(3,1)*m(0,2)*m(1,3) - m(3,1)*m(0,3)*m(1,2);
108  d[3] = -m(0,1)*m(1,2)*m(2,3) + m(0,1)*m(1,3)*m(2,2) + m(1,1)*m(0,2)*m(2,3)
109  - m(1,1)*m(0,3)*m(2,2) - m(2,1)*m(0,2)*m(1,3) + m(2,1)*m(0,3)*m(1,2);
110  d[4] = -m(1,0)*m(2,2)*m(3,3) + m(1,0)*m(2,3)*m(3,2) + m(2,0)*m(1,2)*m(3,3)
111  - m(2,0)*m(1,3)*m(3,2) - m(3,0)*m(1,2)*m(2,3) + m(3,0)*m(1,3)*m(2,2);
112  d[5] = m(0,0)*m(2,2)*m(3,3) - m(0,0)*m(2,3)*m(3,2) - m(2,0)*m(0,2)*m(3,3)
113  + m(2,0)*m(0,3)*m(3,2) + m(3,0)*m(0,2)*m(2,3) - m(3,0)*m(0,3)*m(2,2);
114  d[6] = -m(0,0)*m(1,2)*m(3,3) + m(0,0)*m(1,3)*m(3,2) + m(1,0)*m(0,2)*m(3,3)
115  - m(1,0)*m(0,3)*m(3,2) - m(3,0)*m(0,2)*m(1,3) + m(3,0)*m(0,3)*m(1,2);
116  d[7] = m(0,0)*m(1,2)*m(2,3) - m(0,0)*m(1,3)*m(2,2) - m(1,0)*m(0,2)*m(2,3)
117  + m(1,0)*m(0,3)*m(2,2) + m(2,0)*m(0,2)*m(1,3) - m(2,0)*m(0,3)*m(1,2);
118  d[8] = m(1,0)*m(2,1)*m(3,3) - m(1,0)*m(2,3)*m(3,1) - m(2,0)*m(1,1)*m(3,3)
119  + m(2,0)*m(1,3)*m(3,1) + m(3,0)*m(1,1)*m(2,3) - m(3,0)*m(1,3)*m(2,1);
120  d[9] = -m(0,0)*m(2,1)*m(3,3) + m(0,0)*m(2,3)*m(3,1) + m(2,0)*m(0,1)*m(3,3)
121  - m(2,0)*m(0,3)*m(3,1) - m(3,0)*m(0,1)*m(2,3) + m(3,0)*m(0,3)*m(2,1);
122  d[10]= m(0,0)*m(1,1)*m(3,3) - m(0,0)*m(1,3)*m(3,1) - m(1,0)*m(0,1)*m(3,3)
123  + m(1,0)*m(0,3)*m(3,1) + m(3,0)*m(0,1)*m(1,3) - m(3,0)*m(0,3)*m(1,1);
124  d[11]= -m(0,0)*m(1,1)*m(2,3) + m(0,0)*m(1,3)*m(2,1) + m(1,0)*m(0,1)*m(2,3)
125  - m(1,0)*m(0,3)*m(2,1) - m(2,0)*m(0,1)*m(1,3) + m(2,0)*m(0,3)*m(1,1);
126  d[12]= -m(1,0)*m(2,1)*m(3,2) + m(1,0)*m(2,2)*m(3,1) + m(2,0)*m(1,1)*m(3,2)
127  - m(2,0)*m(1,2)*m(3,1) - m(3,0)*m(1,1)*m(2,2) + m(3,0)*m(1,2)*m(2,1);
128  d[13]= m(0,0)*m(2,1)*m(3,2) - m(0,0)*m(2,2)*m(3,1) - m(2,0)*m(0,1)*m(3,2)
129  + m(2,0)*m(0,2)*m(3,1) + m(3,0)*m(0,1)*m(2,2) - m(3,0)*m(0,2)*m(2,1);
130  d[14]= -m(0,0)*m(1,1)*m(3,2) + m(0,0)*m(1,2)*m(3,1) + m(1,0)*m(0,1)*m(3,2)
131  - m(1,0)*m(0,2)*m(3,1) - m(3,0)*m(0,1)*m(1,2) + m(3,0)*m(0,2)*m(1,1);
132  d[15]= m(0,0)*m(1,1)*m(2,2) - m(0,0)*m(1,2)*m(2,1) - m(1,0)*m(0,1)*m(2,2)
133  + m(1,0)*m(0,2)*m(2,1) + m(2,0)*m(0,1)*m(1,2) - m(2,0)*m(0,2)*m(1,1);
134  return vnl_matrix_fixed<T,4,4>(d);
135 }
136 
137 //: Calculates adjugate_fixed of a small vnl_matrix_fixed (not using svd)
138 // This allows you to write e.g.
139 //
140 // x = vnl_adjugate_fixed(A) * b;
141 //
142 // Note that this function is inlined (except for the call to vnl_det()),
143 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
144 // since that one is using svd.
145 //
146 // \relatesalso vnl_matrix
147 
148 template <class T>
149 vnl_matrix<T> vnl_adjugate_asfixed(vnl_matrix<T> const& m)
150 {
151  assert(m.rows() == m.columns());
152  assert(m.rows() <= 4);
153  if (m.rows() == 1)
154  return vnl_matrix<T>(1,1, T(1)/m(0,0));
155  else if (m.rows() == 2)
156  return vnl_adjugate(vnl_matrix_fixed<T,2,2>(m)).as_ref();
157  else if (m.rows() == 3)
158  return vnl_adjugate(vnl_matrix_fixed<T,3,3>(m)).as_ref();
159  else
160  return vnl_adjugate(vnl_matrix_fixed<T,4,4>(m)).as_ref();
161 }
162 
163 //: Calculates transpose of the adjugate_fixed of a small vnl_matrix_fixed (not using svd)
164 // This allows you to write e.g.
165 //
166 // x = vnl_cofactor(A) * b;
167 //
168 // Note that this function is inlined (except for the call to vnl_det()),
169 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
170 // since that one is using svd. This is also faster than using
171 //
172 // x = vnl_adjugate_fixed(A).transpose() * b;
173 //
174 // \relatesalso vnl_matrix_fixed
175 
176 template <class T>
177 vnl_matrix_fixed<T,1,1> vnl_cofactor(vnl_matrix_fixed<T,1,1> const& m)
178 {
179  return vnl_matrix_fixed<T,1,1>(T(1)/m(0,0));
180 }
181 
182 //: Calculates transpose of the adjugate_fixed of a small vnl_matrix_fixed (not using svd)
183 // This allows you to write e.g.
184 //
185 // x = vnl_cofactor(A) * b;
186 //
187 // Note that this function is inlined (except for the call to vnl_det()),
188 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
189 // since that one is using svd. This is also faster than using
190 //
191 // x = vnl_adjugate_fixed(A).transpose() * b;
192 //
193 // \relatesalso vnl_matrix_fixed
194 
195 template <class T>
196 vnl_matrix_fixed<T,2,2> vnl_cofactor(vnl_matrix_fixed<T,2,2> const& m)
197 {
198 
199  T d[4];
200  d[0] = m(1,1); d[2] = - m(0,1);
201  d[3] = m(0,0); d[1] = - m(1,0);
202  return vnl_matrix_fixed<T,2,2>(d);
203 }
204 
205 //: Calculates transpose of the adjugate_fixed of a small vnl_matrix_fixed (not using svd)
206 // This allows you to write e.g.
207 //
208 // x = vnl_cofactor(A) * b;
209 //
210 // Note that this function is inlined (except for the call to vnl_det()),
211 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
212 // since that one is using svd. This is also faster than using
213 //
214 // x = vnl_adjugate_fixed(A).transpose() * b;
215 //
216 // \relatesalso vnl_matrix_fixed
217 
218 template <class T>
219 vnl_matrix_fixed<T,3,3> vnl_cofactor(vnl_matrix_fixed<T,3,3> const& m)
220 {
221 
222  T d[9];
223  d[0] = (m(1,1)*m(2,2)-m(1,2)*m(2,1));
224  d[3] = (m(2,1)*m(0,2)-m(2,2)*m(0,1));
225  d[6] = (m(0,1)*m(1,2)-m(0,2)*m(1,1));
226  d[1] = (m(1,2)*m(2,0)-m(1,0)*m(2,2));
227  d[4] = (m(0,0)*m(2,2)-m(0,2)*m(2,0));
228  d[7] = (m(1,0)*m(0,2)-m(1,2)*m(0,0));
229  d[2] = (m(1,0)*m(2,1)-m(1,1)*m(2,0));
230  d[5] = (m(0,1)*m(2,0)-m(0,0)*m(2,1));
231  d[8] = (m(0,0)*m(1,1)-m(0,1)*m(1,0));
232  return vnl_matrix_fixed<T,3,3>(d);
233 }
234 
235 //: Calculates transpose of the adjugate_fixed of a small vnl_matrix_fixed (not using svd)
236 // This allows you to write e.g.
237 //
238 // x = vnl_cofactor(A) * b;
239 //
240 // Note that this function is inlined (except for the call to vnl_det()),
241 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
242 // since that one is using svd. This is also faster than using
243 //
244 // x = vnl_adjugate_fixed(A).transpose() * b;
245 //
246 // \relatesalso vnl_matrix_fixed
247 
248 template <class T>
249 vnl_matrix_fixed<T,4,4> vnl_cofactor(vnl_matrix_fixed<T,4,4> const& m)
250 {
251  T d[16];
252  d[0] = m(1,1)*m(2,2)*m(3,3) - m(1,1)*m(2,3)*m(3,2) - m(2,1)*m(1,2)*m(3,3)
253  + m(2,1)*m(1,3)*m(3,2) + m(3,1)*m(1,2)*m(2,3) - m(3,1)*m(1,3)*m(2,2);
254  d[4] = -m(0,1)*m(2,2)*m(3,3) + m(0,1)*m(2,3)*m(3,2) + m(2,1)*m(0,2)*m(3,3)
255  - m(2,1)*m(0,3)*m(3,2) - m(3,1)*m(0,2)*m(2,3) + m(3,1)*m(0,3)*m(2,2);
256  d[8] = m(0,1)*m(1,2)*m(3,3) - m(0,1)*m(1,3)*m(3,2) - m(1,1)*m(0,2)*m(3,3)
257  + m(1,1)*m(0,3)*m(3,2) + m(3,1)*m(0,2)*m(1,3) - m(3,1)*m(0,3)*m(1,2);
258  d[12]= -m(0,1)*m(1,2)*m(2,3) + m(0,1)*m(1,3)*m(2,2) + m(1,1)*m(0,2)*m(2,3)
259  - m(1,1)*m(0,3)*m(2,2) - m(2,1)*m(0,2)*m(1,3) + m(2,1)*m(0,3)*m(1,2);
260  d[1] = -m(1,0)*m(2,2)*m(3,3) + m(1,0)*m(2,3)*m(3,2) + m(2,0)*m(1,2)*m(3,3)
261  - m(2,0)*m(1,3)*m(3,2) - m(3,0)*m(1,2)*m(2,3) + m(3,0)*m(1,3)*m(2,2);
262  d[5] = m(0,0)*m(2,2)*m(3,3) - m(0,0)*m(2,3)*m(3,2) - m(2,0)*m(0,2)*m(3,3)
263  + m(2,0)*m(0,3)*m(3,2) + m(3,0)*m(0,2)*m(2,3) - m(3,0)*m(0,3)*m(2,2);
264  d[9] = -m(0,0)*m(1,2)*m(3,3) + m(0,0)*m(1,3)*m(3,2) + m(1,0)*m(0,2)*m(3,3)
265  - m(1,0)*m(0,3)*m(3,2) - m(3,0)*m(0,2)*m(1,3) + m(3,0)*m(0,3)*m(1,2);
266  d[13]= m(0,0)*m(1,2)*m(2,3) - m(0,0)*m(1,3)*m(2,2) - m(1,0)*m(0,2)*m(2,3)
267  + m(1,0)*m(0,3)*m(2,2) + m(2,0)*m(0,2)*m(1,3) - m(2,0)*m(0,3)*m(1,2);
268  d[2] = m(1,0)*m(2,1)*m(3,3) - m(1,0)*m(2,3)*m(3,1) - m(2,0)*m(1,1)*m(3,3)
269  + m(2,0)*m(1,3)*m(3,1) + m(3,0)*m(1,1)*m(2,3) - m(3,0)*m(1,3)*m(2,1);
270  d[6] = -m(0,0)*m(2,1)*m(3,3) + m(0,0)*m(2,3)*m(3,1) + m(2,0)*m(0,1)*m(3,3)
271  - m(2,0)*m(0,3)*m(3,1) - m(3,0)*m(0,1)*m(2,3) + m(3,0)*m(0,3)*m(2,1);
272  d[10]= m(0,0)*m(1,1)*m(3,3) - m(0,0)*m(1,3)*m(3,1) - m(1,0)*m(0,1)*m(3,3)
273  + m(1,0)*m(0,3)*m(3,1) + m(3,0)*m(0,1)*m(1,3) - m(3,0)*m(0,3)*m(1,1);
274  d[14]= -m(0,0)*m(1,1)*m(2,3) + m(0,0)*m(1,3)*m(2,1) + m(1,0)*m(0,1)*m(2,3)
275  - m(1,0)*m(0,3)*m(2,1) - m(2,0)*m(0,1)*m(1,3) + m(2,0)*m(0,3)*m(1,1);
276  d[3] = -m(1,0)*m(2,1)*m(3,2) + m(1,0)*m(2,2)*m(3,1) + m(2,0)*m(1,1)*m(3,2)
277  - m(2,0)*m(1,2)*m(3,1) - m(3,0)*m(1,1)*m(2,2) + m(3,0)*m(1,2)*m(2,1);
278  d[7] = m(0,0)*m(2,1)*m(3,2) - m(0,0)*m(2,2)*m(3,1) - m(2,0)*m(0,1)*m(3,2)
279  + m(2,0)*m(0,2)*m(3,1) + m(3,0)*m(0,1)*m(2,2) - m(3,0)*m(0,2)*m(2,1);
280  d[11]= -m(0,0)*m(1,1)*m(3,2) + m(0,0)*m(1,2)*m(3,1) + m(1,0)*m(0,1)*m(3,2)
281  - m(1,0)*m(0,2)*m(3,1) - m(3,0)*m(0,1)*m(1,2) + m(3,0)*m(0,2)*m(1,1);
282  d[15]= m(0,0)*m(1,1)*m(2,2) - m(0,0)*m(1,2)*m(2,1) - m(1,0)*m(0,1)*m(2,2)
283  + m(1,0)*m(0,2)*m(2,1) + m(2,0)*m(0,1)*m(1,2) - m(2,0)*m(0,2)*m(1,1);
284  return vnl_matrix_fixed<T,4,4>(d);
285 }
286 
287 //: Calculates transpose of the adjugate_fixed of a small vnl_matrix_fixed (not using svd)
288 // This allows you to write e.g.
289 //
290 // x = vnl_cofactor(A) * b;
291 //
292 // Note that this function is inlined (except for the call to vnl_det()),
293 // which makes it much faster than the vnl_matrix_adjugate_fixed class in vnl/algo
294 // since that one is using svd. This is also faster than using
295 //
296 // x = vnl_adjugate_fixed(A).transpose() * b;
297 //
298 // \relatesalso vnl_matrix
299 
300 template <class T>
301 vnl_matrix<T> vnl_cofactor(vnl_matrix<T> const& m)
302 {
303  assert(m.rows() == m.columns());
304  assert(m.rows() <= 4);
305  if (m.rows() == 1)
306  return vnl_matrix<T>(1,1, T(1)/m(0,0));
307  else if (m.rows() == 2)
308  return vnl_cofactor(vnl_matrix_fixed<T,2,2>(m)).as_ref();
309  else if (m.rows() == 3)
310  return vnl_cofactor(vnl_matrix_fixed<T,3,3>(m)).as_ref();
311  else
312  return vnl_cofactor(vnl_matrix_fixed<T,4,4>(m)).as_ref();
313 }
314 
315 
316 template <class T>
317 vnl_vector_fixed<T,3> vnl_cofactor_row1(vnl_vector_fixed<T,3> const& row2, vnl_vector_fixed<T,3> const& row3)
318 {
319  T d[3];
320  d[0] = (row2[1]*row3[2]-row2[2]*row3[1]);
321  d[1] = (row2[2]*row3[0]-row2[0]*row3[2]);
322  d[2] = (row2[0]*row3[1]-row2[1]*row3[0]);
323  return vnl_vector_fixed<T,3>(d);
324 }
325 
326 #endif // vnl_adjugate_fixed_h_
vnl_vector_fixed< T, 3 > vnl_cofactor_row1(vnl_vector_fixed< T, 3 > const &row2, vnl_vector_fixed< T, 3 > const &row3)
vnl_matrix_fixed< T, 1, 1 > vnl_adjugate(vnl_matrix_fixed< T, 1, 1 > const &m)
vnl_matrix< T > vnl_adjugate_asfixed(vnl_matrix< T > const &m)
vnl_matrix_fixed< T, 1, 1 > vnl_cofactor(vnl_matrix_fixed< T, 1, 1 > const &m)


Generated on 04-09-2015 for elastix by doxygen 1.8.9.1 elastix logo