1 // This is core/vgl/algo/vgl_h_matrix_3d.hxx
2 #ifndef vgl_h_matrix_3d_hxx_
3 #define vgl_h_matrix_3d_hxx_
4
5 #include <iostream>
6 #include <fstream>
7 #include <cmath>
8 #include <limits>
9 #include <cstdlib>
10 #include "vgl_h_matrix_3d.h"
11 #include <cassert>
12 #include <vcl_compiler_detection.h>
13 #ifdef _MSC_VER
14 # include <vcl_msvc_warnings.h>
15 #endif
16 #include <vgl/vgl_plane_3d.h>
17 #include <vnl/vnl_inverse.h>
18 #include <vnl/vnl_numeric_traits.h>
19 #include <vnl/vnl_vector_fixed.h>
20 #include <vnl/vnl_quaternion.h>
21 #include <vnl/algo/vnl_svd.h>
22 # include <vcl_deprecated.h>
23
24 template <class T>
vgl_h_matrix_3d(std::vector<vgl_homg_point_3d<T>> const & points1,std::vector<vgl_homg_point_3d<T>> const & points2)25 vgl_h_matrix_3d<T>::vgl_h_matrix_3d(std::vector<vgl_homg_point_3d<T> > const& points1,
26 std::vector<vgl_homg_point_3d<T> > const& points2)
27 {
28 vnl_matrix<T> W;
29 assert(points1.size() == points2.size());
30 unsigned int numpoints = static_cast<int>( points1.size());
31 if (numpoints < 5)
32 {
33 std::cerr << "\nvhl_h_matrix_3d - minimum of 5 points required\n";
34 std::exit(0);
35 }
36
37 W.set_size(3*numpoints, 16);
38
39 for (unsigned int i = 0; i < numpoints; i++)
40 {
41 T x1 = points1[i].x(), y1 = points1[i].y(), z1 = points1[i].z(), w1 = points1[i].w();
42 T x2 = points2[i].x(), y2 = points2[i].y(), z2 = points2[i].z(), w2 = points2[i].w();
43
44 W[i*3][0]=x1*w2; W[i*3][1]=y1*w2; W[i*3][2]=z1*w2; W[i*3][3]=w1*w2;
45 W[i*3][4]=0.0; W[i*3][5]=0.0; W[i*3][6]=0.0; W[i*3][7]=0.0;
46 W[i*3][8]=0.0; W[i*3][9]=0.0; W[i*3][10]=0.0; W[i*3][11]=0.0;
47 W[i*3][12]=-x1*x2; W[i*3][13]=-y1*x2; W[i*3][14]=-z1*x2; W[i*3][15]=-w1*x2;
48
49 W[i*3+1][0]=0.0; W[i*3+1][1]=0.0; W[i*3+1][2]=0.0; W[i*3+1][3]=0.0;
50 W[i*3+1][4]=x1*w2; W[i*3+1][5]=y1*w2; W[i*3+1][6]=z1*w2; W[i*3+1][7]=w1*w2;
51 W[i*3+1][8]=0.0; W[i*3+1][9]=0.0; W[i*3+1][10]=0.0; W[i*3+1][11]=0.0;
52 W[i*3+1][12]=-x1*y2; W[i*3+1][13]=-y1*y2; W[i*3+1][14]=-z1*y2; W[i*3+1][15]=-w1*y2;
53
54 W[i*3+2][0]=0.0; W[i*3+2][1]=0.0; W[i*3+2][2]=0.0; W[i*3+2][3]=0.0;
55 W[i*3+2][4]=0.0; W[i*3+2][5]=0.0; W[i*3+2][6]=0.0; W[i*3+2][7]=0.0;
56 W[i*3+2][8]=x1*w2; W[i*3+2][9]=y1*w2; W[i*3+2][10]=z1*w2; W[i*3+2][11]=w1*w2;
57 W[i*3+2][12]=-x1*z2; W[i*3+2][13]=-y1*z2; W[i*3+2][14]=-z1*z2; W[i*3+2][15]=-w1*z2;
58 }
59
60 vnl_svd<T> SVD(W);
61 t12_matrix_ = vnl_matrix_fixed<T,4,4>(SVD.nullvector().data_block()); // 16-dim. nullvector
62 }
63
64 template <class T>
vgl_h_matrix_3d(std::istream & s)65 vgl_h_matrix_3d<T>::vgl_h_matrix_3d(std::istream& s)
66 {
67 t12_matrix_.read_ascii(s);
68 }
69
70 template <class T>
vgl_h_matrix_3d(char const * filename)71 vgl_h_matrix_3d<T>::vgl_h_matrix_3d(char const* filename)
72 {
73 std::ifstream f(filename);
74 if (!f.good())
75 std::cerr << "vgl_h_matrix_3d::read: Error opening " << filename << std::endl;
76 else
77 t12_matrix_.read_ascii(f);
78 }
79
80 template <class T>
vgl_h_matrix_3d(vnl_matrix_fixed<T,3,3> const & M,vnl_vector_fixed<T,3> const & m)81 vgl_h_matrix_3d<T>::vgl_h_matrix_3d(vnl_matrix_fixed<T,3,3> const& M,
82 vnl_vector_fixed<T,3> const& m)
83 {
84 for (int r = 0; r < 3; ++r) {
85 for (int c = 0; c < 3; ++c)
86 (t12_matrix_)(r, c) = M(r,c);
87 (t12_matrix_)(r, 3) = m(r);
88 }
89 for (int c = 0; c < 3; ++c)
90 (t12_matrix_)(3,c) = 0;
91 (t12_matrix_)(3,3) = 1;
92 }
93
94 // == OPERATIONS ==
95
96 //-----------------------------------------------------------------------------
97 //
98 template <class T>
99 vgl_homg_point_3d<T>
operator ()(vgl_homg_point_3d<T> const & p) const100 vgl_h_matrix_3d<T>::operator()(vgl_homg_point_3d<T> const& p) const
101 {
102 vnl_vector_fixed<T,4> v2 = t12_matrix_ * vnl_vector_fixed<T,4>(p.x(), p.y(), p.z(), p.w());
103 return vgl_homg_point_3d<T>(v2[0], v2[1], v2[2], v2[3]);
104 }
105
106 template <class T>
operator ()(vgl_pointset_3d<T> const & ptset) const107 vgl_pointset_3d<T> vgl_h_matrix_3d<T>::operator()(vgl_pointset_3d<T> const& ptset) const{
108 vgl_pointset_3d<T> ret;
109 bool has_norms = ptset.has_normals();
110 unsigned np = ptset.npts();
111 std::vector<vgl_point_3d<T> > pts;
112 std::vector<vgl_vector_3d<T> > normals;
113 for(unsigned i =0; i<np; ++i){
114 vgl_homg_point_3d<T> hp = (*this)(vgl_homg_point_3d<T>(ptset.p(i)));
115 pts.push_back(vgl_point_3d<T>(hp));
116 if(has_norms){
117 vgl_vector_3d<T> norm = ptset.n(i);
118 vgl_homg_point_3d<T> hn(norm.x(), norm.y(), norm.z(), 0.0);//direction vector
119 vgl_homg_point_3d<T> hhn = (*this)(hn);
120 vgl_vector_3d<T> hv(hhn.x(), hhn.y(), hhn.z());
121 normals.push_back(hv);
122 }
123 }
124 if(has_norms)
125 ret.set_points_with_normals(pts, normals);
126 else
127 ret.set_points(pts);
128 return ret;
129 }
130
131
132 template <class T>
133 vgl_homg_plane_3d<T>
correlation(vgl_homg_point_3d<T> const & p) const134 vgl_h_matrix_3d<T>::correlation(vgl_homg_point_3d<T> const& p) const
135 {
136 vnl_vector_fixed<T,4> v2 = t12_matrix_ * vnl_vector_fixed<T,4>(p.x(), p.y(), p.z(), p.w());
137 return vgl_homg_plane_3d<T>(v2[0], v2[1], v2[2], v2[3]);
138 }
139
140 template <class T>
141 vgl_homg_plane_3d<T>
preimage(vgl_homg_plane_3d<T> const & l) const142 vgl_h_matrix_3d<T>::preimage(vgl_homg_plane_3d<T> const& l) const
143 {
144 vnl_vector_fixed<T,4> v2 = t12_matrix_.transpose() * vnl_vector_fixed<T,4>(l.a(), l.b(), l.c(), l.d());
145 return vgl_homg_plane_3d<T>(v2[0], v2[1], v2[2], v2[3]);
146 }
147
148 template <class T>
149 vgl_homg_point_3d<T>
correlation(vgl_homg_plane_3d<T> const & l) const150 vgl_h_matrix_3d<T>::correlation(vgl_homg_plane_3d<T> const& l) const
151 {
152 vnl_vector_fixed<T,4> v2 = t12_matrix_ * vnl_vector_fixed<T,4>(l.a(), l.b(), l.c(), l.d());
153 return vgl_homg_point_3d<T>(v2[0], v2[1], v2[2], v2[3]);
154 }
155
156 template <class T>
157 vgl_homg_point_3d<T>
preimage(vgl_homg_point_3d<T> const & p) const158 vgl_h_matrix_3d<T>::preimage(vgl_homg_point_3d<T> const& p) const
159 {
160 vnl_vector_fixed<T,4> v = vnl_inverse(t12_matrix_) * vnl_vector_fixed<T,4>(p.x(), p.y(), p.z(), p.w());
161 return vgl_homg_point_3d<T>(v[0], v[1], v[2], v[3]);
162 }
163 template <class T>
preimage(vgl_pointset_3d<T> const & ptset) const164 vgl_pointset_3d<T> vgl_h_matrix_3d<T>::preimage(vgl_pointset_3d<T> const& ptset) const{
165 vgl_h_matrix_3d<T> hinv = this->get_inverse();
166 return hinv(ptset);
167 }
168
169 template <class T>
170 vgl_homg_plane_3d<T>
operator ()(vgl_homg_plane_3d<T> const & l) const171 vgl_h_matrix_3d<T>::operator()(vgl_homg_plane_3d<T> const& l) const
172 {
173 vnl_vector_fixed<T,4> v = vnl_inverse_transpose(t12_matrix_) * vnl_vector_fixed<T,4>(l.a(), l.b(), l.c(), l.d());
174 return vgl_homg_plane_3d<T>(v[0], v[1], v[2], v[3]);
175 }
176
177 template <class T>
operator <<(std::ostream & s,vgl_h_matrix_3d<T> const & h)178 std::ostream& operator<<(std::ostream& s, vgl_h_matrix_3d<T> const& h)
179 {
180 return s << h.get_matrix();
181 }
182
183 template <class T>
read(std::istream & s)184 bool vgl_h_matrix_3d<T>::read(std::istream& s)
185 {
186 t12_matrix_.read_ascii(s);
187 return s.good() || s.eof();
188 }
189
190 template <class T>
read(char const * filename)191 bool vgl_h_matrix_3d<T>::read(char const* filename)
192 {
193 std::ifstream f(filename);
194 if (!f.good())
195 std::cerr << "vgl_h_matrix_3d::read: Error opening " << filename << std::endl;
196 return read(f);
197 }
198
199 // == DATA ACCESS ==
200
201 template <class T>
get(unsigned int row_index,unsigned int col_index) const202 T vgl_h_matrix_3d<T>::get (unsigned int row_index, unsigned int col_index) const
203 {
204 return t12_matrix_.get(row_index, col_index);
205 }
206
207 template <class T>
get(T * H) const208 void vgl_h_matrix_3d<T>::get (T* H) const
209 {
210 for (T const* iter = t12_matrix_.begin(); iter < t12_matrix_.end(); ++iter)
211 *H++ = *iter;
212 }
213
214 template <class T>
get(vnl_matrix_fixed<T,4,4> * H) const215 void vgl_h_matrix_3d<T>::get (vnl_matrix_fixed<T,4,4>* H) const
216 {
217 *H = t12_matrix_;
218 }
219
220 template <class T>
get(vnl_matrix<T> * H) const221 void vgl_h_matrix_3d<T>::get (vnl_matrix<T>* H) const
222 {
223 VXL_DEPRECATED_MACRO("vgl_h_matrix_3d<T>::get(vnl_matrix<T>*) const");
224 *H = t12_matrix_.as_ref(); // size 4x4
225 }
226
227 template <class T>
228 vgl_h_matrix_3d<T>
get_inverse() const229 vgl_h_matrix_3d<T>::get_inverse() const
230 {
231 return vgl_h_matrix_3d<T>(vnl_inverse(t12_matrix_));
232 }
233
234 template <class T>
235 vgl_h_matrix_3d<T>&
set(T const * H)236 vgl_h_matrix_3d<T>::set (T const* H)
237 {
238 for (T* iter = t12_matrix_.begin(); iter < t12_matrix_.end(); ++iter)
239 *iter = *H++;
240 return *this;
241 }
242
243 template <class T>
244 vgl_h_matrix_3d<T>&
set(vnl_matrix_fixed<T,4,4> const & H)245 vgl_h_matrix_3d<T>::set (vnl_matrix_fixed<T,4,4> const& H)
246 {
247 t12_matrix_ = H;
248 return *this;
249 }
250
251 template <class T>
projective_basis(std::vector<vgl_homg_point_3d<T>> const &)252 bool vgl_h_matrix_3d<T>::projective_basis(std::vector<vgl_homg_point_3d<T> > const& /*five_points*/)
253 {
254 std::cerr << "vgl_h_matrix_3d<T>::projective_basis(5pts) not yet implemented\n";
255 return false;
256 }
257
258 template <class T>
projective_basis(std::vector<vgl_homg_plane_3d<T>> const &)259 bool vgl_h_matrix_3d<T>::projective_basis(std::vector<vgl_homg_plane_3d<T> > const& /*five_planes*/)
260 {
261 std::cerr << "vgl_h_matrix_3d<T>::projective_basis(5planes) not yet implemented\n";
262 return false;
263 }
264
265 template <class T>
266 vgl_h_matrix_3d<T>&
set_identity()267 vgl_h_matrix_3d<T>::set_identity ()
268 {
269 t12_matrix_.set_identity();
270 return *this;
271 }
272
273 template <class T>
274 vgl_h_matrix_3d<T>&
set_translation(T tx,T ty,T tz)275 vgl_h_matrix_3d<T>::set_translation(T tx, T ty, T tz)
276 {
277 t12_matrix_(0, 3) = tx;
278 t12_matrix_(1, 3) = ty;
279 t12_matrix_(2, 3) = tz;
280 return *this;
281 }
282
283 template <class T>
284 vgl_h_matrix_3d<T>&
set_scale(T scale)285 vgl_h_matrix_3d<T>::set_scale(T scale)
286 {
287 for (unsigned r = 0; r<3; ++r)
288 for (unsigned c = 0; c<4; ++c)
289 t12_matrix_[r][c]*=scale;
290 return *this;
291 }
292
293 template <class T>
294 vgl_h_matrix_3d<T>&
set_affine(vnl_matrix_fixed<T,3,4> const & M34)295 vgl_h_matrix_3d<T>::set_affine(vnl_matrix_fixed<T,3,4> const& M34)
296 {
297 for (unsigned r = 0; r<3; ++r)
298 for (unsigned c = 0; c<4; ++c)
299 t12_matrix_[r][c] = M34[r][c];
300 t12_matrix_[3][0] = t12_matrix_[3][1] = t12_matrix_[3][2] = T(0); t12_matrix_[3][3] = T(1);
301 return *this;
302 }
303
304 template <class T>
305 vgl_h_matrix_3d<T>&
set_rotation_about_axis(vnl_vector_fixed<T,3> const & axis,T angle)306 vgl_h_matrix_3d<T>::set_rotation_about_axis(vnl_vector_fixed<T,3> const& axis, T angle)
307 {
308 vnl_quaternion<T> q(axis, angle);
309 //get the transpose of the rotation matrix
310 vnl_matrix_fixed<T,3,3> R = q.rotation_matrix_transpose();
311 //fill in with the transpose
312 for (int c = 0; c<3; c++)
313 for (int r = 0; r<3; r++)
314 t12_matrix_[r][c]=R[c][r];
315 return *this;
316 }
317
318 template <class T>
319 vgl_h_matrix_3d<T>&
set_rotation_roll_pitch_yaw(T yaw,T pitch,T roll)320 vgl_h_matrix_3d<T>::set_rotation_roll_pitch_yaw(T yaw, T pitch, T roll)
321 {
322 typedef typename vnl_numeric_traits<T>::real_t real_t;
323 real_t ax = yaw/2, ay = pitch/2, az = roll/2;
324
325 vnl_quaternion<T> qx((T)std::sin(ax),0,0,(T)std::cos(ax));
326 vnl_quaternion<T> qy(0,(T)std::sin(ay),0,(T)std::cos(ay));
327 vnl_quaternion<T> qz(0,0,(T)std::sin(az),(T)std::cos(az));
328 vnl_quaternion<T> q = qz*qy*qx;
329
330 vnl_matrix_fixed<T,3,3> R = q.rotation_matrix_transpose();
331 //fill in with the transpose
332 for (int c = 0; c<3; c++)
333 for (int r = 0; r<3; r++)
334 t12_matrix_[r][c]=R[c][r];
335 return *this;
336 }
337
338 template <class T>
339 vgl_h_matrix_3d<T>&
set_rotation_euler(T rz1,T ry,T rz2)340 vgl_h_matrix_3d<T>::set_rotation_euler(T rz1, T ry, T rz2)
341 {
342 typedef typename vnl_numeric_traits<T>::real_t real_t;
343 real_t az1 = rz1/2, ay = ry/2, az2 = rz2/2;
344
345 vnl_quaternion<T> qz1(0,0,T(std::sin(az1)),T(std::cos(az1)));
346 vnl_quaternion<T> qy(0,T(std::sin(ay)),0,T(std::cos(ay)));
347 vnl_quaternion<T> qz2(0,0,T(std::sin(az2)),T(std::cos(az2)));
348 vnl_quaternion<T> q = qz2*qy*qz1;
349
350 vnl_matrix_fixed<T,3,3> R = q.rotation_matrix_transpose();
351 //fill in with the transpose
352 for (int c = 0; c<3; c++)
353 for (int r = 0; r<3; r++)
354 t12_matrix_[r][c]=R[c][r];
355 return *this;
356 }
357
358 template <class T>
359 vgl_h_matrix_3d<T>&
set_rotation_matrix(vnl_matrix_fixed<T,3,3> const & R)360 vgl_h_matrix_3d<T>::set_rotation_matrix(vnl_matrix_fixed<T,3,3> const& R)
361 {
362 for (unsigned r = 0; r<3; ++r)
363 for (unsigned c = 0; c<3; ++c)
364 t12_matrix_[r][c] = R[r][c];
365 return *this;
366 }
367
368
369 template <class T>
370 void
set_reflection_plane(vgl_plane_3d<double> const & l)371 vgl_h_matrix_3d<T>::set_reflection_plane(vgl_plane_3d<double> const& l)
372 {
373 t12_matrix_.fill(T(0));
374 t12_matrix_(0,0) = T(l.nx()*l.nx());
375 t12_matrix_(1,1) = T(l.ny()*l.ny());
376 t12_matrix_(2,2) = T(l.nz()*l.nz());
377 t12_matrix_(0,1) = t12_matrix_(1,0) = T(l.nx()*l.ny());
378 t12_matrix_(0,2) = t12_matrix_(2,0) = T(l.nx()*l.nz());
379 t12_matrix_(1,2) = t12_matrix_(2,1) = T(l.ny()*l.nz());
380 t12_matrix_(0,3) = T(l.nx()*l.d());
381 t12_matrix_(1,3) = T(l.ny()*l.d());
382 t12_matrix_(2,3) = T(l.nz()*l.d());
383 t12_matrix_ *= -2/(t12_matrix_(0,0)+t12_matrix_(1,1)+t12_matrix_(2,2));
384 t12_matrix_(0,0) += (T)1;
385 t12_matrix_(1,1) += (T)1;
386 t12_matrix_(2,2) += (T)1;
387 t12_matrix_(3,3) += (T)1;
388 }
389
390
391 template <class T>
is_rotation() const392 bool vgl_h_matrix_3d<T>::is_rotation() const
393 {
394 return t12_matrix_.get(0,3) == (T)0
395 && t12_matrix_.get(1,3) == (T)0
396 && t12_matrix_.get(2,3) == (T)0
397 && this->is_euclidean();
398 }
399
400
401 template <class T>
is_euclidean() const402 bool vgl_h_matrix_3d<T>::is_euclidean() const
403 {
404 T eps = 10*std::numeric_limits<T>::epsilon();
405 if ( t12_matrix_.get(3,0) != (T)0 ||
406 t12_matrix_.get(3,1) != (T)0 ||
407 t12_matrix_.get(3,2) != (T)0 ||
408 std::fabs(t12_matrix_.get(3,3)-T(1)) > eps)
409 return false; // should not have a projective part
410
411 // use an error tolerance on the orthonormality constraint
412 vnl_matrix_fixed<T,3,3> R = get_upper_3x3_matrix();
413 R *= R.transpose();
414 R(0,0) -= T(1);
415 R(1,1) -= T(1);
416 R(2,2) -= T(1);
417 return R.absolute_value_max() <= eps;
418 }
419
420 template <class T>
is_affine() const421 bool vgl_h_matrix_3d<T>::is_affine() const{
422 if ( t12_matrix_.get(3,0) != (T)0 ||
423 t12_matrix_.get(3,1) != (T)0 ||
424 t12_matrix_.get(3,2) != (T)0 ||
425 std::fabs(t12_matrix_.get(3,3)) > 10*std::numeric_limits<T>::epsilon())
426 return false; // should not have a projective part
427 return !(this->is_euclidean());
428 }
429
430 template <class T>
is_identity() const431 bool vgl_h_matrix_3d<T>::is_identity() const
432 {
433 return t12_matrix_.is_identity();
434 }
435
436
437 template <class T>
438 vgl_h_matrix_3d<T>
get_upper_3x3() const439 vgl_h_matrix_3d<T>::get_upper_3x3() const
440 {
441 //only sensible for affine transformations
442 T d = t12_matrix_[3][3];
443 assert(d<-1e-9 || d>1e-9);
444 vnl_matrix_fixed<T,4,4> m(0.0);
445 for (unsigned r = 0; r<3; r++)
446 for (unsigned c = 0; c<3; c++)
447 m[r][c] = t12_matrix_[r][c]/d;
448 m[3][3]=1.0;
449 return vgl_h_matrix_3d<T>(m);
450 }
451
452 template <class T>
453 vnl_matrix_fixed<T,3,3>
get_upper_3x3_matrix() const454 vgl_h_matrix_3d<T>::get_upper_3x3_matrix() const
455 {
456 vnl_matrix_fixed<T,3,3> R;
457 vgl_h_matrix_3d<T> m = this->get_upper_3x3();
458 for (unsigned r = 0; r<3; r++)
459 for (unsigned c = 0; c<3; c++)
460 R[r][c] = m.get(r,c);
461 return R;
462 }
463 template <class T>
polar_decomposition(vnl_matrix_fixed<T,3,3> & S,vnl_matrix_fixed<T,3,3> & R) const464 void vgl_h_matrix_3d<T>::polar_decomposition(vnl_matrix_fixed<T, 3, 3>& S, vnl_matrix_fixed<T, 3, 3>& R) const{
465 vnl_matrix_fixed<T, 3, 3> up = this->get_upper_3x3_matrix();
466 vnl_matrix<T> M{up.as_matrix()};
467 vnl_svd<T> svd(M);
468 vnl_matrix<T> U = svd.U();
469 vnl_matrix<T> W{ svd.W().as_matrix() };
470 vnl_matrix<T> V = svd.V();
471 R = vnl_matrix_fixed<T, 3, 3> ( U*V.transpose());
472 S = vnl_matrix_fixed<T, 3, 3> (V*W*V.transpose());
473 return;
474 }
475
476 template <class T>
477 vgl_homg_point_3d<T>
get_translation() const478 vgl_h_matrix_3d<T>::get_translation() const
479 {
480 //only sensible for affine transformations
481 T d = t12_matrix_[3][3];
482 assert(d<-1e-9 || d>1e-9);
483 return vgl_homg_point_3d<T>(t12_matrix_[0][3]/d,
484 t12_matrix_[1][3]/d,
485 t12_matrix_[2][3]/d,
486 (T)1);
487 }
488
489 template <class T>
490 vnl_vector_fixed<T,3>
get_translation_vector() const491 vgl_h_matrix_3d<T>::get_translation_vector() const
492 {
493 vgl_homg_point_3d<T> p = this->get_translation();
494 return vnl_vector_fixed<T,3>(p.x(), p.y(), p.z());
495 }
496
497 //----------------------------------------------------------------------------
498 #undef VGL_H_MATRIX_3D_INSTANTIATE
499 #define VGL_H_MATRIX_3D_INSTANTIATE(T) \
500 template class vgl_h_matrix_3d<T >; \
501 template std::ostream& operator<<(std::ostream&, vgl_h_matrix_3d<T > const& ); \
502 template std::istream& operator>>(std::istream&, vgl_h_matrix_3d<T >& )
503
504 #endif // vgl_h_matrix_3d_hxx_
505