Point Cloud Library (PCL) 1.14.0
Loading...
Searching...
No Matches
auxiliary.h
1/*
2 * Software License Agreement (BSD License)
3 *
4 * Point Cloud Library (PCL) - www.pointclouds.org
5 * Copyright (c) 2012-, Open Perception, Inc.
6 *
7 * All rights reserved.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 *
13 * * Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * * Redistributions in binary form must reproduce the above
16 * copyright notice, this list of conditions and the following
17 * disclaimer in the documentation and/or other materials provided
18 * with the distribution.
19 * * Neither the name of the copyright holder(s) nor the names of its
20 * contributors may be used to endorse or promote products derived
21 * from this software without specific prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
33 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 * POSSIBILITY OF SUCH DAMAGE.
35 *
36 */
37
38#pragma once
39
40#include <Eigen/Core> // for Matrix
41#include <Eigen/Geometry> // for AngleAxis
42#include <pcl/point_types.h>
43
44#define AUX_PI_FLOAT 3.14159265358979323846f
45#define AUX_HALF_PI 1.57079632679489661923f
46#define AUX_DEG_TO_RADIANS (3.14159265358979323846f/180.0f)
47
48namespace pcl
49{
50 namespace recognition
51 {
52 namespace aux
53 {
54 template<typename T> bool
55 compareOrderedPairs (const std::pair<T,T>& a, const std::pair<T,T>& b)
56 {
57 if ( a.first == b.first )
58 return a.second < b.second;
59
60 return a.first < b.first;
61 }
62
63 template<typename T> T
64 sqr (T a)
65 {
66 return (a*a);
67 }
68
69 template<typename T> T
70 clamp (T value, T min, T max)
71 {
72 if ( value < min )
73 return min;
74 if ( value > max )
75 return max;
76
77 return value;
78 }
79
80 /** \brief Expands the destination bounding box 'dst' such that it contains 'src'. */
81 template<typename T> void
82 expandBoundingBox (T dst[6], const T src[6])
83 {
84 if ( src[0] < dst[0] ) dst[0] = src[0];
85 if ( src[2] < dst[2] ) dst[2] = src[2];
86 if ( src[4] < dst[4] ) dst[4] = src[4];
87
88 if ( src[1] > dst[1] ) dst[1] = src[1];
89 if ( src[3] > dst[3] ) dst[3] = src[3];
90 if ( src[5] > dst[5] ) dst[5] = src[5];
91 }
92
93 /** \brief Expands the bounding box 'bbox' such that it contains the point 'p'. */
94 template<typename T> void
95 expandBoundingBoxToContainPoint (T bbox[6], const T p[3])
96 {
97 if ( p[0] < bbox[0] ) bbox[0] = p[0];
98 else if ( p[0] > bbox[1] ) bbox[1] = p[0];
99
100 if ( p[1] < bbox[2] ) bbox[2] = p[1];
101 else if ( p[1] > bbox[3] ) bbox[3] = p[1];
102
103 if ( p[2] < bbox[4] ) bbox[4] = p[2];
104 else if ( p[2] > bbox[5] ) bbox[5] = p[2];
105 }
106
107 /** \brief v[0] = v[1] = v[2] = value */
108 template <typename T> void
109 set3 (T v[3], T value)
110 {
111 v[0] = v[1] = v[2] = value;
112 }
113
114 /** \brief dst = src */
115 template <typename T> void
116 copy3 (const T src[3], T dst[3])
117 {
118 dst[0] = src[0];
119 dst[1] = src[1];
120 dst[2] = src[2];
121 }
122
123 /** \brief dst = src */
124 template <typename T> void
125 copy3 (const T src[3], pcl::PointXYZ& dst)
126 {
127 dst.x = src[0];
128 dst.y = src[1];
129 dst.z = src[2];
130 }
131
132 /** \brief a = -a */
133 template <typename T> void
134 flip3 (T a[3])
135 {
136 a[0] = -a[0];
137 a[1] = -a[1];
138 a[2] = -a[2];
139 }
140
141 /** \brief a = b */
142 template <typename T> bool
143 equal3 (const T a[3], const T b[3])
144 {
145 return (a[0] == b[0] && a[1] == b[1] && a[2] == b[2]);
146 }
147
148 /** \brief a += b */
149 template <typename T> void
150 add3 (T a[3], const T b[3])
151 {
152 a[0] += b[0];
153 a[1] += b[1];
154 a[2] += b[2];
155 }
156
157 /** \brief c = a + b */
158 template <typename T> void
159 sum3 (const T a[3], const T b[3], T c[3])
160 {
161 c[0] = a[0] + b[0];
162 c[1] = a[1] + b[1];
163 c[2] = a[2] + b[2];
164 }
165
166 /** \brief c = a - b */
167 template <typename T> void
168 diff3 (const T a[3], const T b[3], T c[3])
169 {
170 c[0] = a[0] - b[0];
171 c[1] = a[1] - b[1];
172 c[2] = a[2] - b[2];
173 }
174
175 template <typename T> void
176 cross3 (const T v1[3], const T v2[3], T out[3])
177 {
178 out[0] = v1[1]*v2[2] - v1[2]*v2[1];
179 out[1] = v1[2]*v2[0] - v1[0]*v2[2];
180 out[2] = v1[0]*v2[1] - v1[1]*v2[0];
181 }
182
183 /** \brief Returns the length of v. */
184 template <typename T> T
185 length3 (const T v[3])
186 {
187 return (std::sqrt (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]));
188 }
189
190 /** \brief Returns the Euclidean distance between a and b. */
191 template <typename T> T
192 distance3 (const T a[3], const T b[3])
193 {
194 T l[3] = {a[0]-b[0], a[1]-b[1], a[2]-b[2]};
195 return (length3 (l));
196 }
197
198 /** \brief Returns the squared Euclidean distance between a and b. */
199 template <typename T> T
200 sqrDistance3 (const T a[3], const T b[3])
201 {
202 return (aux::sqr (a[0]-b[0]) + aux::sqr (a[1]-b[1]) + aux::sqr (a[2]-b[2]));
203 }
204
205 /** \brief Returns the dot product a*b */
206 template <typename T> T
207 dot3 (const T a[3], const T b[3])
208 {
209 return (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]);
210 }
211
212 /** \brief Returns the dot product (x, y, z)*(u, v, w) = x*u + y*v + z*w */
213 template <typename T> T
214 dot3 (T x, T y, T z, T u, T v, T w)
215 {
216 return (x*u + y*v + z*w);
217 }
218
219 /** \brief v = scalar*v. */
220 template <typename T> void
221 mult3 (T* v, T scalar)
222 {
223 v[0] *= scalar;
224 v[1] *= scalar;
225 v[2] *= scalar;
226 }
227
228 /** \brief out = scalar*v. */
229 template <typename T> void
230 mult3 (const T* v, T scalar, T* out)
231 {
232 out[0] = v[0]*scalar;
233 out[1] = v[1]*scalar;
234 out[2] = v[2]*scalar;
235 }
236
237 /** \brief Normalize v */
238 template <typename T> void
239 normalize3 (T v[3])
240 {
241 T inv_len = (static_cast<T> (1.0))/aux::length3 (v);
242 v[0] *= inv_len;
243 v[1] *= inv_len;
244 v[2] *= inv_len;
245 }
246
247 /** \brief Returns the square length of v. */
248 template <typename T> T
249 sqrLength3 (const T v[3])
250 {
251 return (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
252 }
253
254 /** Projects 'x' on the plane through 0 and with normal 'planeNormal' and saves the result in 'out'. */
255 template <typename T> void
256 projectOnPlane3 (const T x[3], const T planeNormal[3], T out[3])
257 {
258 T dot = aux::dot3 (planeNormal, x);
259 // Project 'x' on the plane normal
260 T nproj[3] = {-dot*planeNormal[0], -dot*planeNormal[1], -dot*planeNormal[2]};
261 aux::sum3 (x, nproj, out);
262 }
263
264 /** \brief Sets 'm' to the 3x3 identity matrix. */
265 template <typename T> void
266 identity3x3 (T m[9])
267 {
268 m[0] = m[4] = m[8] = 1.0;
269 m[1] = m[2] = m[3] = m[5] = m[6] = m[7] = 0.0;
270 }
271
272 /** \brief out = mat*v. 'm' is an 1D array of 9 elements treated as a 3x3 matrix (row major order). */
273 template <typename T> void
274 mult3x3(const T m[9], const T v[3], T out[3])
275 {
276 out[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
277 out[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
278 out[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
279 }
280
281 /** Let x, y, z be the columns of the matrix a = [x|y|z]. The method computes out = a*m.
282 * Note that 'out' is a 1D array of 9 elements and the resulting matrix is stored in row
283 * major order, i.e., the first matrix row is (out[0] out[1] out[2]), the second
284 * (out[3] out[4] out[5]) and the third (out[6] out[7] out[8]). */
285 template <typename T> void
286 mult3x3 (const T x[3], const T y[3], const T z[3], const T m[3][3], T out[9])
287 {
288 out[0] = x[0]*m[0][0] + y[0]*m[1][0] + z[0]*m[2][0];
289 out[1] = x[0]*m[0][1] + y[0]*m[1][1] + z[0]*m[2][1];
290 out[2] = x[0]*m[0][2] + y[0]*m[1][2] + z[0]*m[2][2];
291
292 out[3] = x[1]*m[0][0] + y[1]*m[1][0] + z[1]*m[2][0];
293 out[4] = x[1]*m[0][1] + y[1]*m[1][1] + z[1]*m[2][1];
294 out[5] = x[1]*m[0][2] + y[1]*m[1][2] + z[1]*m[2][2];
295
296 out[6] = x[2]*m[0][0] + y[2]*m[1][0] + z[2]*m[2][0];
297 out[7] = x[2]*m[0][1] + y[2]*m[1][1] + z[2]*m[2][1];
298 out[8] = x[2]*m[0][2] + y[2]*m[1][2] + z[2]*m[2][2];
299 }
300
301 /** \brief The first 9 elements of 't' are treated as a 3x3 matrix (row major order) and the last 3 as a translation.
302 * First, 'p' is multiplied by that matrix and then translated. The result is saved in 'out'. */
303 template<class T> void
304 transform(const T t[12], const T p[3], T out[3])
305 {
306 out[0] = t[0]*p[0] + t[1]*p[1] + t[2]*p[2] + t[9];
307 out[1] = t[3]*p[0] + t[4]*p[1] + t[5]*p[2] + t[10];
308 out[2] = t[6]*p[0] + t[7]*p[1] + t[8]*p[2] + t[11];
309 }
310
311 /** \brief The first 9 elements of 't' are treated as a 3x3 matrix (row major order) and the last 3 as a translation.
312 * First, (x, y, z) is multiplied by that matrix and then translated. The result is saved in 'out'. */
313 template<class T> void
314 transform(const T t[12], T x, T y, T z, T out[3])
315 {
316 out[0] = t[0]*x + t[1]*y + t[2]*z + t[9];
317 out[1] = t[3]*x + t[4]*y + t[5]*z + t[10];
318 out[2] = t[6]*x + t[7]*y + t[8]*z + t[11];
319 }
320
321 /** \brief Compute out = (upper left 3x3 of mat)*p + last column of mat. */
322 template<class T> void
323 transform(const Eigen::Matrix<T,4,4>& mat, const pcl::PointXYZ& p, pcl::PointXYZ& out)
324 {
325 out.x = mat(0,0)*p.x + mat(0,1)*p.y + mat(0,2)*p.z + mat(0,3);
326 out.y = mat(1,0)*p.x + mat(1,1)*p.y + mat(1,2)*p.z + mat(1,3);
327 out.z = mat(2,0)*p.x + mat(2,1)*p.y + mat(2,2)*p.z + mat(2,3);
328 }
329
330 /** \brief The first 9 elements of 't' are treated as a 3x3 matrix (row major order) and the last 3 as a translation.
331 * First, 'p' is multiplied by that matrix and then translated. The result is saved in 'out'. */
332 template<class T> void
333 transform(const T t[12], const pcl::PointXYZ& p, T out[3])
334 {
335 out[0] = t[0]*p.x + t[1]*p.y + t[2]*p.z + t[9];
336 out[1] = t[3]*p.x + t[4]*p.y + t[5]*p.z + t[10];
337 out[2] = t[6]*p.x + t[7]*p.y + t[8]*p.z + t[11];
338 }
339
340 /** \brief Returns true if the points 'p1' and 'p2' are co-planar and false otherwise. The method assumes that 'n1'
341 * is a normal at 'p1' and 'n2' is a normal at 'p2'. 'max_angle' is the threshold used for the test. The bigger
342 * the value the larger the deviation between the normals can be which still leads to a positive test result. The
343 * angle has to be in radians. */
344 template<typename T> bool
345 pointsAreCoplanar (const T p1[3], const T n1[3], const T p2[3], const T n2[3], T max_angle)
346 {
347 // Compute the angle between 'n1' and 'n2' and compare it with 'max_angle'
348 if ( std::acos (aux::clamp (aux::dot3 (n1, n2), -1.0f, 1.0f)) > max_angle )
349 return (false);
350
351 T cl[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
352 aux::normalize3 (cl);
353
354 // Compute the angle between 'cl' and 'n1'
355 T tmp_angle = std::acos (aux::clamp (aux::dot3 (n1, cl), -1.0f, 1.0f));
356
357 // 'tmp_angle' should not deviate too much from 90 degrees
358 if ( std::fabs (tmp_angle - AUX_HALF_PI) > max_angle )
359 return (false);
360
361 // All tests passed => the points are coplanar
362 return (true);
363 }
364
365 template<typename Scalar> void
366 array12ToMatrix4x4 (const Scalar src[12], Eigen::Matrix<Scalar, 4, 4>& dst)
367 {
368 dst(0,0) = src[0]; dst(0,1) = src[1]; dst(0,2) = src[2]; dst(0,3) = src[9];
369 dst(1,0) = src[3]; dst(1,1) = src[4]; dst(1,2) = src[5]; dst(1,3) = src[10];
370 dst(2,0) = src[6]; dst(2,1) = src[7]; dst(2,2) = src[8]; dst(2,3) = src[11];
371 dst(3,0) = dst(3,1) = dst(3,2) = 0.0; dst(3,3) = 1.0;
372 }
373
374 template<typename Scalar> void
375 matrix4x4ToArray12 (const Eigen::Matrix<Scalar, 4, 4>& src, Scalar dst[12])
376 {
377 dst[0] = src(0,0); dst[1] = src(0,1); dst[2] = src(0,2); dst[9] = src(0,3);
378 dst[3] = src(1,0); dst[4] = src(1,1); dst[5] = src(1,2); dst[10] = src(1,3);
379 dst[6] = src(2,0); dst[7] = src(2,1); dst[8] = src(2,2); dst[11] = src(2,3);
380 }
381
382 /** \brief The method copies the input array 'src' to the eigen matrix 'dst' in row major order.
383 * dst[0] = src(0,0); dst[1] = src(0,1); dst[2] = src(0,2);
384 * dst[3] = src(1,0); dst[4] = src(1,1); dst[5] = src(1,2);
385 * dst[6] = src(2,0); dst[7] = src(2,1); dst[8] = src(2,2);
386 * */
387 template <typename T> void
388 eigenMatrix3x3ToArray9RowMajor (const Eigen::Matrix<T,3,3>& src, T dst[9])
389 {
390 dst[0] = src(0,0); dst[1] = src(0,1); dst[2] = src(0,2);
391 dst[3] = src(1,0); dst[4] = src(1,1); dst[5] = src(1,2);
392 dst[6] = src(2,0); dst[7] = src(2,1); dst[8] = src(2,2);
393 }
394
395 /** \brief The method copies the input array 'src' to the eigen matrix 'dst' in row major order.
396 * dst(0,0) = src[0]; dst(0,1) = src[1]; dst(0,2) = src[2];
397 * dst(1,0) = src[3]; dst(1,1) = src[4]; dst(1,2) = src[5];
398 * dst(2,0) = src[6]; dst(2,1) = src[7]; dst(2,2) = src[8];
399 * */
400 template <typename T> void
401 toEigenMatrix3x3RowMajor (const T src[9], Eigen::Matrix<T,3,3>& dst)
402 {
403 dst(0,0) = src[0]; dst(0,1) = src[1]; dst(0,2) = src[2];
404 dst(1,0) = src[3]; dst(1,1) = src[4]; dst(1,2) = src[5];
405 dst(2,0) = src[6]; dst(2,1) = src[7]; dst(2,2) = src[8];
406 }
407
408 /** brief Computes a rotation matrix from the provided input vector 'axis_angle'. The direction of 'axis_angle' is the rotation axis
409 * and its magnitude is the angle of rotation about that axis. 'rotation_matrix' is the output rotation matrix saved in row major order. */
410 template <typename T> void
411 axisAngleToRotationMatrix (const T axis_angle[3], T rotation_matrix[9])
412 {
413 // Get the angle of rotation
414 T angle = aux::length3 (axis_angle);
415 if ( angle == 0.0 )
416 {
417 // Undefined rotation -> set to identity
418 aux::identity3x3 (rotation_matrix);
419 return;
420 }
421
422 // Normalize the input
423 T normalized_axis_angle[3];
424 aux::mult3 (axis_angle, static_cast<T> (1.0)/angle, normalized_axis_angle);
425
426 // The eigen objects
427 Eigen::Matrix<T,3,1> mat_axis(normalized_axis_angle);
428 Eigen::AngleAxis<T> eigen_angle_axis (angle, mat_axis);
429
430 // Save the output
431 aux::eigenMatrix3x3ToArray9RowMajor (eigen_angle_axis.toRotationMatrix (), rotation_matrix);
432 }
433
434 /** brief Extracts the angle-axis representation from 'rotation_matrix', i.e., computes a rotation 'axis' and an 'angle'
435 * of rotation about that axis from the provided input. The output 'angle' is in the range [0, pi] and 'axis' is normalized. */
436 template <typename T> void
437 rotationMatrixToAxisAngle (const T rotation_matrix[9], T axis[3], T& angle)
438 {
439 // The eigen objects
440 Eigen::AngleAxis<T> angle_axis;
441 Eigen::Matrix<T,3,3> rot_mat;
442 // Copy the input matrix to the eigen matrix in row major order
443 aux::toEigenMatrix3x3RowMajor (rotation_matrix, rot_mat);
444
445 // Do the computation
446 angle_axis.fromRotationMatrix (rot_mat);
447
448 // Save the result
449 axis[0] = angle_axis.axis () (0,0);
450 axis[1] = angle_axis.axis () (1,0);
451 axis[2] = angle_axis.axis () (2,0);
452 angle = angle_axis.angle ();
453
454 // Make sure that 'angle' is in the range [0, pi]
455 if ( angle > AUX_PI_FLOAT )
456 {
457 angle = 2.0f*AUX_PI_FLOAT - angle;
458 aux::flip3 (axis);
459 }
460 }
461 } // namespace aux
462 } // namespace recognition
463} // namespace pcl
Defines all the PCL implemented PointT point type structures.
T sqrLength3(const T v[3])
Returns the square length of v.
Definition auxiliary.h:249
void expandBoundingBoxToContainPoint(T bbox[6], const T p[3])
Expands the bounding box 'bbox' such that it contains the point 'p'.
Definition auxiliary.h:95
void array12ToMatrix4x4(const Scalar src[12], Eigen::Matrix< Scalar, 4, 4 > &dst)
Definition auxiliary.h:366
T distance3(const T a[3], const T b[3])
Returns the Euclidean distance between a and b.
Definition auxiliary.h:192
void axisAngleToRotationMatrix(const T axis_angle[3], T rotation_matrix[9])
brief Computes a rotation matrix from the provided input vector 'axis_angle'.
Definition auxiliary.h:411
T dot3(const T a[3], const T b[3])
Returns the dot product a*b.
Definition auxiliary.h:207
T sqrDistance3(const T a[3], const T b[3])
Returns the squared Euclidean distance between a and b.
Definition auxiliary.h:200
bool equal3(const T a[3], const T b[3])
a = b
Definition auxiliary.h:143
T length3(const T v[3])
Returns the length of v.
Definition auxiliary.h:185
void mult3(T *v, T scalar)
v = scalar*v.
Definition auxiliary.h:221
void matrix4x4ToArray12(const Eigen::Matrix< Scalar, 4, 4 > &src, Scalar dst[12])
Definition auxiliary.h:375
void add3(T a[3], const T b[3])
a += b
Definition auxiliary.h:150
void transform(const T t[12], const T p[3], T out[3])
The first 9 elements of 't' are treated as a 3x3 matrix (row major order) and the last 3 as a transla...
Definition auxiliary.h:304
void normalize3(T v[3])
Normalize v.
Definition auxiliary.h:239
T clamp(T value, T min, T max)
Definition auxiliary.h:70
void eigenMatrix3x3ToArray9RowMajor(const Eigen::Matrix< T, 3, 3 > &src, T dst[9])
The method copies the input array 'src' to the eigen matrix 'dst' in row major order.
Definition auxiliary.h:388
void toEigenMatrix3x3RowMajor(const T src[9], Eigen::Matrix< T, 3, 3 > &dst)
The method copies the input array 'src' to the eigen matrix 'dst' in row major order.
Definition auxiliary.h:401
void diff3(const T a[3], const T b[3], T c[3])
c = a - b
Definition auxiliary.h:168
void set3(T v[3], T value)
v[0] = v[1] = v[2] = value
Definition auxiliary.h:109
void projectOnPlane3(const T x[3], const T planeNormal[3], T out[3])
Projects 'x' on the plane through 0 and with normal 'planeNormal' and saves the result in 'out'.
Definition auxiliary.h:256
void flip3(T a[3])
a = -a
Definition auxiliary.h:134
bool pointsAreCoplanar(const T p1[3], const T n1[3], const T p2[3], const T n2[3], T max_angle)
Returns true if the points 'p1' and 'p2' are co-planar and false otherwise.
Definition auxiliary.h:345
void copy3(const T src[3], T dst[3])
dst = src
Definition auxiliary.h:116
void sum3(const T a[3], const T b[3], T c[3])
c = a + b
Definition auxiliary.h:159
void cross3(const T v1[3], const T v2[3], T out[3])
Definition auxiliary.h:176
void expandBoundingBox(T dst[6], const T src[6])
Expands the destination bounding box 'dst' such that it contains 'src'.
Definition auxiliary.h:82
void identity3x3(T m[9])
Sets 'm' to the 3x3 identity matrix.
Definition auxiliary.h:266
bool compareOrderedPairs(const std::pair< T, T > &a, const std::pair< T, T > &b)
Definition auxiliary.h:55
void mult3x3(const T m[9], const T v[3], T out[3])
out = mat*v.
Definition auxiliary.h:274
void rotationMatrixToAxisAngle(const T rotation_matrix[9], T axis[3], T &angle)
brief Extracts the angle-axis representation from 'rotation_matrix', i.e., computes a rotation 'axis'...
Definition auxiliary.h:437
A point structure representing Euclidean xyz coordinates.