Loading...
Searching...
No Matches
analyticalGeometry.h
Go to the documentation of this file.
// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
#ifndef BALL_MATHS_ANALYTICALGEOMETRY_H
#define BALL_MATHS_ANALYTICALGEOMETRY_H
#ifndef BALL_COMMON_EXCEPTION_H
# include <BALL/COMMON/exception.h>
#endif
#ifndef BALL_MATHS_ANGLE_H
# include <BALL/MATHS/angle.h>
#endif
#ifndef BALL_MATHS_CIRCLE3_H
# include <BALL/MATHS/circle3.h>
#endif
#ifndef BALL_MATHS_LINE3_H
# include <BALL/MATHS/line3.h>
#endif
#ifndef BALL_MATHS_PLANE3_H
# include <BALL/MATHS/plane3.h>
#endif
#ifndef BALL_MATHS_SPHERE3_H
# include <BALL/MATHS/sphere3.h>
#endif
#ifndef BALL_MATHS_VECTOR3_H
# include <BALL/MATHS/vector3.h>
#endif
#define BALL_MATRIX_CELL(m, dim, row, col) *((m) + (row) * (dim) + (col))
#define BALL_CELL(x, y) *((m) + (y) * (dim) + (x))
namespace BALL
{
template <typename T>
BALL_INLINE
T getDeterminant_(const T* m, Size dim)
{
T determinant = 0;
Index dim1 = dim - 1;
if (dim > 1)
{
T* submatrix = new T[dim1 * dim1];
for (Index i = 0; i < (Index)dim; ++i)
{
for (Index j = 0; j < dim1; ++j)
{
for (Index k = 0; k < dim1; ++k)
{
*(submatrix + j * dim1 + k) = *(m + (j + 1) * dim + (k < i ? k : k + 1));
}
}
determinant += *(m + i) * (i / 2.0 == i / 2 ? 1 : -1) * getDeterminant_(submatrix, dim1);
}
delete [] submatrix;
}
else
{
determinant = *m;
}
return determinant;
}
template <typename T>
T getDeterminant(const T* m, Size dim)
{
if (dim == 2)
{
return (BALL_CELL(0,0) * BALL_CELL(1,1) - BALL_CELL(0,1) * BALL_CELL(1,0));
}
else if (dim == 3)
{
return ( BALL_CELL(0,0) * BALL_CELL(1,1) * BALL_CELL(2,2)
+ BALL_CELL(0,1) * BALL_CELL(1,2) * BALL_CELL(2,0)
+ BALL_CELL(0,2) * BALL_CELL(1,0) * BALL_CELL(2,1)
- BALL_CELL(0,2) * BALL_CELL(1,1) * BALL_CELL(2,0)
- BALL_CELL(0,0) * BALL_CELL(1,2) * BALL_CELL(2,1)
- BALL_CELL(0,1) * BALL_CELL(1,0) * BALL_CELL(2,2));
}
else
{
return getDeterminant_(m, dim);
}
}
template <typename T>
BALL_INLINE
T getDeterminant2(const T* m)
{
Size dim = 2;
return (BALL_CELL(0,0) * BALL_CELL(1,1) - BALL_CELL(0,1) * BALL_CELL(1,0));
}
template <typename T>
BALL_INLINE
T getDeterminant2(const T& m00, const T& m01, const T& m10, const T& m11)
{
return (m00 * m11 - m01 * m10);
}
template <typename T>
BALL_INLINE
T getDeterminant3(const T *m)
{
Size dim = 3;
return ( BALL_CELL(0,0) * BALL_CELL(1,1) * BALL_CELL(2,2)
+ BALL_CELL(0,1) * BALL_CELL(1,2) * BALL_CELL(2,0)
+ BALL_CELL(0,2) * BALL_CELL(1,0) * BALL_CELL(2,1)
- BALL_CELL(0,2) * BALL_CELL(1,1) * BALL_CELL(2,0)
- BALL_CELL(0,0) * BALL_CELL(1,2) * BALL_CELL(2,1)
- BALL_CELL(0,1) * BALL_CELL(1,0) * BALL_CELL(2,2));
}
template <typename T>
BALL_INLINE T
getDeterminant3(const T& m00, const T& m01, const T& m02, const T& m10, const T& m11, const T& m12, const T& m20, const T& m21, const T& m22)
{
return ( m00 * m11 * m22 + m01 * m12 * m20 + m02 * m10 * m21 - m02 * m11 * m20 - m00 * m12 * m21 - m01 * m10 * m22);
}
template <typename T>
bool SolveSystem(const T* m, T* x, const Size dim)
{
T pivot;
Index i, j, k, p;
// the column dimension of the matrix
const Size col_dim = dim + 1;
T* matrix = new T[dim * (dim + 1)];
const T* source = m;
T* target = (T*)matrix;
T* end = (T*)&BALL_MATRIX_CELL(matrix, col_dim, dim - 1, dim);
while (target <= end)
{
*target++ = *source++;
}
for (i = 0; i < (Index)dim; ++i)
{
pivot = BALL_MATRIX_CELL(matrix, col_dim, i, i);
p = i;
for (j = i + 1; j < (Index)dim; ++j)
{
if (Maths::isLess(pivot, BALL_MATRIX_CELL(matrix, col_dim, j, i)))
{
pivot = BALL_MATRIX_CELL(matrix, col_dim, j, i);
p = j;
}
}
if (p != i)
{
T tmp;
for (k = i; k < (Index)dim + 1; ++k)
{
tmp = BALL_MATRIX_CELL(matrix, dim, i, k);
BALL_MATRIX_CELL(matrix, col_dim, i, k) = BALL_MATRIX_CELL(matrix, col_dim, p, k);
BALL_MATRIX_CELL(matrix, col_dim, p, k) = tmp;
}
}
else if (Maths::isZero(pivot) || Maths::isNan(pivot))
{
// invariant: matrix m is singular
delete [] matrix;
return false;
}
for (j = dim; j >= i; --j)
{
BALL_MATRIX_CELL(matrix, col_dim, i, j) /= pivot;
}
for (j = i + 1; j < (Index)dim; ++j)
{
pivot = BALL_MATRIX_CELL(matrix, col_dim, j, i);
for (k = dim; k>= i; --k)
{
BALL_MATRIX_CELL(matrix, col_dim, j, k) -= pivot * BALL_MATRIX_CELL(matrix, col_dim, i, k);
}
}
}
x[dim - 1] = BALL_MATRIX_CELL(matrix, col_dim, dim - 1, dim);
for (i = dim - 2; i >= 0; --i)
{
x[i] = BALL_MATRIX_CELL(matrix, col_dim, i, dim);
for (j = i + 1; j < (Index)dim; ++j)
{
x[i] -= BALL_MATRIX_CELL(matrix, col_dim, i, j) * x[j];
}
}
delete [] matrix;
return true;
}
#undef BALL_CELL
#undef BALL_MATRIX_CELL
template <typename T>
BALL_INLINE
bool SolveSystem2(const T& a1, const T& b1, const T& c1, const T& a2, const T& b2, const T& c2, T& x1, T& x2)
{
T quot = (a1 * b2 - a2 * b1);
if (Maths::isZero(quot))
{
return false;
}
x1 = (c1 * b2 - c2 * b1) / quot;
x2 = (a1 * c2 - a2 * c1) / quot;
return true;
}
template <typename T>
short SolveQuadraticEquation(const T& a, const T& b, const T &c, T &x1, T &x2)
{
if (a == 0)
{
if (b == 0)
{
return 0;
}
x1 = x2 = c / b;
return 1;
}
T discriminant = b * b - 4 * a * c;
if (Maths::isLess(discriminant, 0))
{
return 0;
}
T sqrt_discriminant = sqrt(discriminant);
if (Maths::isZero(sqrt_discriminant))
{
x1 = x2 = -b / (2 * a);
return 1;
}
else
{
x1 = (-b + sqrt_discriminant) / (2 * a);
x2 = (-b - sqrt_discriminant) / (2 * a);
return 2;
}
}
template <typename T>
BALL_INLINE
TVector3<T> GetPartition(const TVector3<T>& a, const TVector3<T>& b)
{
return TVector3<T>((b.x + a.x) / 2, (b.y + a.y) / 2, (b.z + a.z) / 2);
}
template <typename T>
BALL_INLINE
TVector3<T> GetPartition(const TVector3<T>& a, const TVector3<T>& b, const T& r, const T& s)
{
T sum = r + s;
if (sum == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return TVector3<T>
((s * a.x + r * b.x) / sum,
(s * a.y + r * b.y) / sum,
(s * a.z + r * b.z) / sum);
}
template <typename T>
BALL_INLINE
T GetDistance(const TVector3<T>& a, const TVector3<T>& b)
{
T dx = a.x - b.x;
T dy = a.y - b.y;
T dz = a.z - b.z;
return sqrt(dx * dx + dy * dy + dz * dz);
}
template <typename T>
BALL_INLINE
T GetDistance(const TLine3<T>& line, const TVector3<T>& point)
{
if (line.d.getLength() == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return ((line.d % (point - line.p)).getLength() / line.d.getLength());
}
template <typename T>
BALL_INLINE
T GetDistance(const TVector3<T>& point, const TLine3<T>& line)
{
return GetDistance(line, point);
}
template <typename T>
T GetDistance(const TLine3<T>& a, const TLine3<T>& b)
{
T cross_product_length = (a.d % b.d).getLength();
if (Maths::isZero(cross_product_length))
{ // invariant: parallel lines
if (a.d.getLength() == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return ((a.d % (b.p - a.p)).getLength() / a.d.getLength());
}
else
{
T spat_product = TVector3<T>::getTripleProduct(a.d, b.d, b.p - a.p);
if (Maths::isNotZero(spat_product))
{ // invariant: windschiefe lines
return (Maths::abs(spat_product) / cross_product_length);
}
else
{
// invariant: intersecting lines
return 0;
}
}
}
template <typename T>
BALL_INLINE
T GetDistance(const TVector3<T>& point, const TPlane3<T>& plane)
{
T length = plane.n.getLength();
if (length == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return (Maths::abs(plane.n * (point - plane.p)) / length);
}
template <typename T>
BALL_INLINE
T GetDistance(const TPlane3<T>& plane, const TVector3<T>& point)
{
return GetDistance(point, plane);
}
template <typename T>
BALL_INLINE
T GetDistance(const TLine3<T>& line, const TPlane3<T>& plane)
{
T length = plane.n.getLength();
if (length == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return (Maths::abs(plane.n * (line.p - plane.p)) / length);
}
template <typename T>
BALL_INLINE
T GetDistance(const TPlane3<T>& plane, const TLine3<T>& line)
{
return GetDistance(line, plane);
}
template <typename T>
BALL_INLINE
T GetDistance(const TPlane3<T>& a, const TPlane3<T>& b)
{
T length = a.n.getLength();
if (length == (T)0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
return (Maths::abs(a.n * (a.p - b.p)) / length);
}
template <typename T>
BALL_INLINE
bool GetAngle(const TVector3<T>& a, const TVector3<T>& b, TAngle<T> &intersection_angle)
{
T length_product = a.getSquareLength() * b.getSquareLength();
if(Maths::isZero(length_product))
{
return false;
}
intersection_angle = a.getAngle(b);
return true;
}
template <typename T>
BALL_INLINE
bool GetAngle(const TLine3<T>& a, const TLine3<T>& b, TAngle<T>& intersection_angle)
{
T length_product = a.d.getSquareLength() * b.d.getSquareLength();
if(Maths::isZero(length_product))
{
return false;
}
intersection_angle = acos(Maths::abs(a.d * b.d) / sqrt(length_product));
return true;
}
template <typename T>
BALL_INLINE
bool GetAngle(const TPlane3<T>& plane, const TVector3<T>& vector, TAngle<T>& intersection_angle)
{
T length_product = plane.n.getSquareLength() * vector.getSquareLength();
if (Maths::isZero(length_product))
{
return false;
}
else
{
intersection_angle = asin(Maths::abs(plane.n * vector) / sqrt(length_product));
return true;
}
}
template <typename T>
BALL_INLINE
bool GetAngle(const TVector3<T>& vector ,const TPlane3<T>& plane, TAngle<T> &intersection_angle)
{
return GetAngle(plane, vector, intersection_angle);
}
template <typename T>
BALL_INLINE
bool GetAngle(const TPlane3<T>& plane,const TLine3<T>& line, TAngle<T>& intersection_angle)
{
T length_product = plane.n.getSquareLength() * line.d.getSquareLength();
if (Maths::isZero(length_product))
{
return false;
}
intersection_angle = asin(Maths::abs(plane.n * line.d) / sqrt(length_product));
return true;
}
template <typename T>
BALL_INLINE
bool GetAngle(const TLine3<T>& line, const TPlane3<T>& plane, TAngle<T>& intersection_angle)
{
return GetAngle(plane, line, intersection_angle);
}
template <typename T>
BALL_INLINE
bool GetAngle(const TPlane3<T>& a, const TPlane3<T>& b, TAngle<T>& intersection_angle)
{
T length_product = a.n.getSquareLength() * b.n.getSquareLength();
if(Maths::isZero(length_product))
{
return false;
}
intersection_angle = acos(Maths::abs(a.n * b.n) / sqrt(length_product));
return true;
}
template <typename T>
bool GetIntersection(const TLine3<T>& a, const TLine3<T>& b, TVector3<T>& point)
{
T c1, c2;
if ((SolveSystem2(a.d.x, -b.d.x, b.p.x - a.p.x, a.d.y, -b.d.y, b.p.y - a.p.y, c1, c2) == true && Maths::isEqual(a.p.z + a.d.z * c1, b.p.z + b.d.z * c2)) || (SolveSystem2(a.d.x, -b.d.x, b.p.x - a.p.x, a.d.z, -b.d.z, b.p.z - a.p.z, c1, c2) == true && Maths::isEqual(a.p.y + a.d.y * c1, b.p.y + b.d.y * c2)) || (SolveSystem2(a.d.y, -b.d.y, b.p.y - a.p.y, a.d.z, -b.d.z, b.p.z - a.p.z, c1, c2) == true && Maths::isEqual(a.p.x + a.d.x * c1, b.p.x + b.d.x * c2)))
{
point.set(a.p.x + a.d.x * c1, a.p.y + a.d.y * c1, a.p.z + a.d.z * c1);
return true;
}
return false;
}
template <typename T>
BALL_INLINE
bool GetIntersection(const TPlane3<T>& plane, const TLine3<T>& line, TVector3<T>& intersection_point)
{
T dot_product = plane.n * line.d;
if (Maths::isZero(dot_product))
{
return false;
}
intersection_point.set(line.p + (plane.n * (plane.p - line.p)) * line.d / dot_product);
return true;
}
template <typename T>
BALL_INLINE
bool GetIntersection(const TLine3<T>& line, const TPlane3<T>& plane, TVector3<T>& intersection_point)
{
return GetIntersection(plane, line, intersection_point);
}
template <typename T>
bool GetIntersection(const TPlane3<T>& plane1, const TPlane3<T>& plane2, TLine3<T>& line)
{
T u = plane1.p*plane1.n;
T v = plane2.p*plane2.n;
T det = plane1.n.x*plane2.n.y-plane1.n.y*plane2.n.x;
if (Maths::isZero(det))
{
det = plane1.n.x*plane2.n.z-plane1.n.z*plane2.n.x;
if (Maths::isZero(det))
{
det = plane1.n.y*plane2.n.z-plane1.n.z*plane2.n.y;
if (Maths::isZero(det))
{
return false;
}
else
{
T a = plane2.n.z/det;
T b = -plane1.n.z/det;
T c = -plane2.n.y/det;
T d = plane1.n.y/det;
line.p.x = 0;
line.p.y = a*u+b*v;
line.p.z = c*u+d*v;
line.d.x = -1;
line.d.y = a*plane1.n.x+b*plane2.n.x;
line.d.z = c*plane1.n.x+d*plane2.n.x;
}
}
else
{
T a = plane2.n.z/det;
T b = -plane1.n.z/det;
T c = -plane2.n.x/det;
T d = plane1.n.x/det;
line.p.x = a*u+b*v;
line.p.y = 0;
line.p.z = c*u+d*v;
line.d.x = a*plane1.n.y+b*plane2.n.y;
line.d.y = -1;
line.d.z = c*plane1.n.y+d*plane2.n.y;
}
}
else
{
T a = plane2.n.y/det;
T b = -plane1.n.y/det;
T c = -plane2.n.x/det;
T d = plane1.n.x/det;
line.p.x = a*u+b*v;
line.p.y = c*u+d*v;
line.p.z = 0;
line.d.x = a*plane1.n.z+b*plane2.n.z;
line.d.y = c*plane1.n.z+d*plane2.n.z;
line.d.z = -1;
}
return true;
}
template <typename T>
bool GetIntersection(const TSphere3<T>& sphere, const TLine3<T>& line, TVector3<T>& intersection_point1, TVector3<T>& intersection_point2)
{
T x1, x2;
short number_of_solutions = SolveQuadraticEquation (line.d * line.d, (line.p - sphere.p) * line.d * 2, (line.p - sphere.p) * (line.p - sphere.p) - sphere.radius * sphere.radius, x1, x2);
if (number_of_solutions == 0)
{
return false;
}
intersection_point1 = line.p + x1 * line.d;
intersection_point2 = line.p + x2 * line.d;
return true;
}
template <typename T>
BALL_INLINE
bool GetIntersection(const TLine3<T>& line, const TSphere3<T>& sphere, TVector3<T>& intersection_point1, TVector3<T>& intersection_point2)
{
return GetIntersection(sphere, line, intersection_point1, intersection_point2);
}
template <typename T>
bool GetIntersection(const TSphere3<T>& sphere, const TPlane3<T>& plane, TCircle3<T>& intersection_circle)
{
T distance = GetDistance(sphere.p, plane);
if (Maths::isGreater(distance, sphere.radius))
{
return false;
}
TVector3<T> Vector3(plane.n);
Vector3.normalize();
if (Maths::isEqual(distance, sphere.radius))
{
intersection_circle.set(sphere.p + sphere.radius * Vector3, plane.n, 0);
}
else
{
intersection_circle.set
(sphere.p + distance * Vector3, plane.n,
sqrt(sphere.radius * sphere.radius - distance * distance));
}
return true;
}
template <typename T>
BALL_INLINE bool
GetIntersection(const TPlane3<T>& plane, const TSphere3<T>& sphere, TCircle3<T>& intersection_circle)
{
return GetIntersection(sphere, plane, intersection_circle);
}
template <typename T>
bool GetIntersection(const TSphere3<T>& a, const TSphere3<T>& b, TCircle3<T>& intersection_circle)
{
TVector3<T> norm = b.p - a.p;
T square_dist = norm * norm;
if (Maths::isZero(square_dist))
{
return false;
}
T dist = sqrt(square_dist);
if (Maths::isLess(a.radius + b.radius, dist))
{
return false;
}
if (Maths::isGreaterOrEqual(Maths::abs(a.radius - b.radius), dist))
{
return false;
}
T radius1_square = a.radius * a.radius;
T radius2_square = b.radius * b.radius;
T u = radius1_square - radius2_square + square_dist;
T length = u / (2 * square_dist);
T square_radius = radius1_square - u * length / 2;
if (square_radius < 0)
{
return false;
}
intersection_circle.p = a.p + (norm * length);
intersection_circle.radius = sqrt(square_radius);
intersection_circle.n = norm / dist;
return true;
}
template <class T>
bool GetIntersection(const TSphere3<T>& s1, const TSphere3<T>& s2, const TSphere3<T>& s3, TVector3<T>& p1, TVector3<T>& p2, bool test = true)
{
T r1_square = s1.radius*s1.radius;
T r2_square = s2.radius*s2.radius;
T r3_square = s3.radius*s3.radius;
T p1_square_length = s1.p*s1.p;
T p2_square_length = s2.p*s2.p;
T p3_square_length = s3.p*s3.p;
T u = (r2_square-r1_square-p2_square_length+p1_square_length)/2;
T v = (r3_square-r1_square-p3_square_length+p1_square_length)/2;
TPlane3<T> plane1;
TPlane3<T> plane2;
try
{
plane1 = TPlane3<T>(s2.p.x-s1.p.x,s2.p.y-s1.p.y,s2.p.z-s1.p.z,u);
plane2 = TPlane3<T>(s3.p.x-s1.p.x,s3.p.y-s1.p.y,s3.p.z-s1.p.z,v);
}
catch (Exception::DivisionByZero&)
{
return false;
}
TLine3<T> line;
if (GetIntersection(plane1,plane2,line))
{
TVector3<T> diff(s1.p-line.p);
T x1 = 0;
T x2 = 0;
if (SolveQuadraticEquation(line.d*line.d, -diff*line.d*2, diff*diff-r1_square, x1,x2) > 0)
{
p1 = line.p+x1*line.d;
p2 = line.p+x2*line.d;
if (test)
{
TVector3<T> test = s1.p-p1;
if (Maths::isNotEqual(test*test,r1_square))
{
return false;
}
test = s1.p-p2;
if (Maths::isNotEqual(test*test,r1_square))
{
return false;
}
test = s2.p-p1;
if (Maths::isNotEqual(test*test,r2_square))
{
return false;
}
test = s2.p-p2;
if (Maths::isNotEqual(test*test,r2_square))
{
return false;
}
test = s3.p-p1;
if (Maths::isNotEqual(test*test,r3_square))
{
return false;
}
test = s3.p-p2;
if (Maths::isNotEqual(test*test,r3_square))
{
return false;
}
}
return true;
}
}
return false;
}
template <typename T>
BALL_INLINE
bool isCollinear(const TVector3<T>& a, const TVector3<T>& b)
{
return (a % b).isZero();
}
template <typename T>
BALL_INLINE
bool isComplanar(const TVector3<T>& a, const TVector3<T>& b, const TVector3<T>& c)
{
return Maths::isZero(TVector3<T>::getTripleProduct(a, b, c));
}
template <typename T>
BALL_INLINE
bool isComplanar(const TVector3<T>& a, const TVector3<T>& b, const TVector3<T>& c, const TVector3<T>& d)
{
return isComplanar(a - b, a - c, a - d);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TVector3<T>& a, const TVector3<T>& b)
{
return Maths::isZero(a * b);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TVector3<T>& vector, const TLine3<T>& line)
{
return Maths::isZero(vector * line.d);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TLine3<T>& line, const TVector3<T>& vector)
{
return isOrthogonal(vector, line);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TLine3<T>& a, const TLine3<T>& b)
{
return Maths::isZero(a.d * b.d);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TVector3<T>& vector, const TPlane3<T>& plane)
{
return isCollinear(vector, plane.n);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TPlane3<T>& plane, const TVector3<T>& vector)
{
return isOrthogonal(vector, plane);
}
template <typename T>
BALL_INLINE
bool isOrthogonal(const TPlane3<T>& a, const TPlane3<T>& b)
{
return Maths::isZero(a.n * b.n);
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TVector3<T>& point, const TLine3<T>& line)
{
return Maths::isZero(GetDistance(point, line));
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TLine3<T>& line, const TVector3<T>& point)
{
return isIntersecting(point, line);
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TLine3<T>& a, const TLine3<T>& b)
{
return Maths::isZero(GetDistance(a, b));
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TVector3<T>& point, const TPlane3<T>& plane)
{
return Maths::isZero(GetDistance(point, plane));
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TPlane3<T>& plane, const TVector3<T>& point)
{
return isIntersecting(point, plane);
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TLine3<T>& line, const TPlane3<T>& plane)
{
return Maths::isZero(GetDistance(line, plane));
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TPlane3<T>& plane, const TLine3<T>& line)
{
return isIntersecting(line, plane);
}
template <typename T>
BALL_INLINE
bool isIntersecting(const TPlane3<T>& a, const TPlane3<T>& b)
{
return Maths::isZero(GetDistance(a, b));
}
template <typename T>
BALL_INLINE
bool isParallel(const TLine3<T>& line, const TPlane3<T>& plane)
{
return isOrthogonal(line.d, plane.n);
}
template <typename T>
BALL_INLINE
bool isParallel(const TPlane3<T>& plane, const TLine3<T>& line)
{
return isParallel(line, plane);
}
template <typename T>
BALL_INLINE
bool isParallel(const TPlane3<T>& a, const TPlane3<T>& b)
{
return isCollinear(a.n, b.n);
}
template <typename T>
TAngle<T> getOrientedAngle
(const T& ax, const T& ay, const T& az,
const T& bx, const T& by, const T& bz,
const T& nx, const T& ny, const T& nz)
{
// Calculate the length of the two normals
T bl = (T) sqrt((double)ax * ax + ay * ay + az * az);
T el = (T) sqrt((double)bx * bx + by * by + bz * bz);
T bel = (T) (ax * bx + ay * by + az * bz);
// if one or both planes are degenerated
if (bl * el == 0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
bel /= (bl * el);
if (bel > 1.0)
{
bel = 1;
}
else if (bel < -1.0)
{
bel = -1;
}
T acosbel = (T) acos(bel); // >= 0
if (( nx * (az * by - ay * bz)
+ ny * (ax * bz - az * bx)
+ nz * (ay * bx - ax * by)) > 0)
{
acosbel = Constants::PI+Constants::PI-acosbel;
}
return TAngle<T>(acosbel);
}
template <typename T>
BALL_INLINE
TAngle<T>getOrientedAngle(const TVector3<T>& a, const TVector3<T>& b, const TVector3<T>& normal)
{
return getOrientedAngle(a.x, a.y, a.z, b.x, b.y, b.z, normal.x, normal.y, normal.z);
}
template <typename T>
TAngle<T> getTorsionAngle
(const T& ax, const T& ay, const T& az,
const T& bx, const T& by, const T& bz,
const T& cx, const T& cy, const T& cz,
const T& dx, const T& dy, const T& dz)
{
T abx = ax - bx;
T aby = ay - by;
T abz = az - bz;
T cbx = cx - bx;
T cby = cy - by;
T cbz = cz - bz;
T cdx = cx - dx;
T cdy = cy - dy;
T cdz = cz - dz;
// Calculate the normals to the two planes n1 and n2
// this is given as the cross products:
// AB x BC
// --------- = n1
// |AB x BC|
//
// BC x CD
// --------- = n2
// |BC x CD|
// Normal to plane 1
T ndax = aby * cbz - abz * cby;
T nday = abz * cbx - abx * cbz;
T ndaz = abx * cby - aby * cbx;
// Normal to plane 2
T neax = cbz * cdy - cby * cdz;
T neay = cbx * cdz - cbz * cdx;
T neaz = cby * cdx - cbx * cdy;
// Calculate the length of the two normals
T bl = (T) sqrt((double)ndax * ndax + nday * nday + ndaz * ndaz);
T el = (T) sqrt((double)neax * neax + neay * neay + neaz * neaz);
T bel = (T) (ndax * neax + nday * neay + ndaz * neaz);
// if one or both planes are degenerated
if (bl * el == 0)
{
throw Exception::DivisionByZero(__FILE__, __LINE__);
}
bel /= (bl * el);
if (bel > 1.0)
{
bel = 1;
}
else if (bel < -1.0)
{
bel = -1;
}
T acosbel = (T) acos(bel);
if ((cbx * (ndaz * neay - nday * neaz)
+ cby * (ndax * neaz - ndaz * neax)
+ cbz * (nday * neax - ndax * neay))
< 0)
{
acosbel = -acosbel;
}
acosbel = (acosbel > 0.0)
? Constants::PI - acosbel
: -(Constants::PI + acosbel);
return TAngle<T>(acosbel);
}
} // namespace BALL
#endif // BALL_MATHS_ANALYTICALGEOMETRY_H
exception.h
BALL_MATRIX_CELL
#define BALL_MATRIX_CELL(m, dim, row, col)
Definition analyticalGeometry.h:36
BALL_CELL
#define BALL_CELL(x, y)
Definition analyticalGeometry.h:37
angle.h
circle3.h
BALL::Exception::DivisionByZero
Definition COMMON/exception.h:342
BALL::TAngle
Definition angle.h:53
BALL::TCircle3
Definition circle3.h:46
BALL::TCircle3::p
TVector3< T > p
Definition circle3.h:277
BALL::TCircle3::radius
T radius
Definition circle3.h:287
BALL::TCircle3::n
TVector3< T > n
Definition circle3.h:282
BALL::TCircle3::set
void set(const TCircle3 &circle)
Definition circle3.h:134
BALL::TLine3
Definition line3.h:47
BALL::TLine3::p
TVector3< T > p
Definition line3.h:347
BALL::TLine3::d
TVector3< T > d
Definition line3.h:351
BALL::TPlane3
Definition plane3.h:53
BALL::TPlane3::p
TVector3< T > p
Definition plane3.h:366
BALL::TPlane3::n
TVector3< T > n
Definition plane3.h:370
BALL::TSphere3
Definition sphere3.h:33
BALL::TSphere3::p
TVector3< T > p
Definition sphere3.h:250
BALL::TSphere3::radius
T radius
Definition sphere3.h:254
BALL::TVector3
Definition vector3.h:71
BALL::TVector3::getAngle
TAngle< T > getAngle(const TVector3 &vector) const
Definition vector3.h:971
BALL::TVector3::x
T x
Definition vector3.h:489
BALL::TVector3::getSquareLength
T getSquareLength() const
Definition vector3.h:764
BALL::TVector3::getTripleProduct
static T getTripleProduct(const TVector3< T > &a, const TVector3< T > &b, const TVector3< T > &c)
Definition vector3.h:1016
BALL::TVector3::y
T y
Definition vector3.h:493
BALL::TVector3::set
void set(const T *ptr)
Definition vector3.h:615
BALL::TVector3::z
T z
Definition vector3.h:497
BALL_INLINE
#define BALL_INLINE
Definition debug.h:15
line3.h
BALL::Constants::PI
BALL_EXTERN_VARIABLE const double PI
PI.
Definition constants.h:35
BALL::Maths::isNan
bool isNan(const T &t)
Definition MATHS/common.h:160
BALL::Maths::isNotZero
bool isNotZero(const T &t)
Definition MATHS/common.h:200
BALL::Maths::isNotEqual
bool isNotEqual(const T1 &a, const T2 &b)
Definition MATHS/common.h:224
BALL::Maths::isGreater
bool isGreater(const T1 &a, const T2 &b)
Definition MATHS/common.h:273
BALL::Maths::isGreaterOrEqual
bool isGreaterOrEqual(const T1 &a, const T2 &b)
Definition MATHS/common.h:261
BALL::Maths::isLess
bool isLess(const T1 &a, const T2 &b)
Definition MATHS/common.h:236
BALL::Maths::isZero
bool isZero(const T &t)
Definition MATHS/common.h:189
BALL::Maths::abs
T abs(const T &t)
Definition MATHS/common.h:50
BALL::Maths::isEqual
bool isEqual(const T1 &a, const T2 &b)
Definition MATHS/common.h:212
BALL
Definition constants.h:13
BALL::SolveSystem2
BALL_INLINE bool SolveSystem2(const T &a1, const T &b1, const T &c1, const T &a2, const T &b2, const T &c2, T &x1, T &x2)
Definition analyticalGeometry.h:285
BALL::getDeterminant_
BALL_INLINE T getDeterminant_(const T *m, Size dim)
Definition analyticalGeometry.h:57
BALL::GetAngle
BALL_INLINE bool GetAngle(const TVector3< T > &a, const TVector3< T > &b, TAngle< T > &intersection_angle)
Definition analyticalGeometry.h:551
BALL::isCollinear
BALL_INLINE bool isCollinear(const TVector3< T > &a, const TVector3< T > &b)
Definition analyticalGeometry.h:1007
BALL::Size
BALL_SIZE_TYPE Size
Definition COMMON/global.h:114
BALL::GetDistance
BALL_INLINE T GetDistance(const TVector3< T > &a, const TVector3< T > &b)
Definition analyticalGeometry.h:387
BALL::SolveQuadraticEquation
short SolveQuadraticEquation(const T &a, const T &b, const T &c, T &x1, T &x2)
Definition analyticalGeometry.h:310
BALL::isOrthogonal
BALL_INLINE bool isOrthogonal(const TVector3< T > &a, const TVector3< T > &b)
Definition analyticalGeometry.h:1046
BALL::GetPartition
BALL_INLINE TVector3< T > GetPartition(const TVector3< T > &a, const TVector3< T > &b)
Definition analyticalGeometry.h:351
BALL::getOrientedAngle
TAngle< T > getOrientedAngle(const T &ax, const T &ay, const T &az, const T &bx, const T &by, const T &bz, const T &nx, const T &ny, const T &nz)
Definition analyticalGeometry.h:1260
BALL::GetIntersection
bool GetIntersection(const TLine3< T > &a, const TLine3< T > &b, TVector3< T > &point)
Definition analyticalGeometry.h:681
BALL::getDeterminant
T getDeterminant(const T *m, Size dim)
Definition analyticalGeometry.h:93
BALL::Vector3
TVector3< float > Vector3
Definition vector3.h:1084
BALL::getTorsionAngle
TAngle< T > getTorsionAngle(const T &ax, const T &ay, const T &az, const T &bx, const T &by, const T &bz, const T &cx, const T &cy, const T &cz, const T &dx, const T &dy, const T &dz)
Definition analyticalGeometry.h:1324
BALL::getDeterminant3
BALL_INLINE T getDeterminant3(const T *m)
Definition analyticalGeometry.h:143
BALL::isComplanar
BALL_INLINE bool isComplanar(const TVector3< T > &a, const TVector3< T > &b, const TVector3< T > &c)
Definition analyticalGeometry.h:1020
BALL::getDeterminant2
BALL_INLINE T getDeterminant2(const T *m)
Definition analyticalGeometry.h:119
BALL::isIntersecting
BALL_INLINE bool isIntersecting(const TVector3< T > &point, const TLine3< T > &line)
Definition analyticalGeometry.h:1130
BALL::SolveSystem
bool SolveSystem(const T *m, T *x, const Size dim)
Definition analyticalGeometry.h:191
BALL::isParallel
BALL_INLINE bool isParallel(const TLine3< T > &line, const TPlane3< T > &plane)
Definition analyticalGeometry.h:1226
BALL::Index
BALL_INDEX_TYPE Index
Definition COMMON/global.h:105
plane3.h
sphere3.h
vector3.h