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