mathlib.h File Reference

#include <cassert>
#include <iostream>
#include <string>
#include <vector>
#include <point.h>
#include <print.h>
#include <typedefs.h>
#include <zero.h>

Include dependency graph for mathlib.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

class crossproduct

In 3D given two vectors that make a plane find a third vector at right angles to this plane. More...

class intervalintersection

Interval overlap tests. More...

class solverInconsistent< T >

Solve 1D and 2D linear equations with one move variable than the number of equations. More...

class surfaceplane< T >

class convexity< T >

class orthoproj< T, V, IP, Dim >

class determinant

Determinant calculations. More...

class solver< T >

class transrotate2D

class circleLine

Circle and line intersection test. More...

Defines

#define PI 3.1415926535897932384626433

Functions

template<typename T >

void unitbound (T &val)

Mape the variable to be inside [0.0,1.0].

void integersetdiff (vector< uint > &vecnot, vector< uint > &vec, uintc N)

Given the integer set 0,1,.

void mathcombination (double &result, uintc n, uintc k)

n!/((n-k)!k!)

void tetrahedronvolume (double &vol, point3< double > const &p0, point3< double > const &p1, point3< double > const &p2, point3< double > const &p3)

Calculate the volume of the tetrahedron.

void polygonconvexarea (double &area, vector< point2< double > > const &v)

Calculate the area of the convex polygon.

void trianglearea (double &f, doublec x0, doublec y0, doublec x1, doublec y1, doublec x2, doublec y2)

void trianglearea (double &f, point3< double > const &p0, point3< double > const &p1, point3< double > const &p2)

Calcluate the triangles area given its three points.

template<typename T >

boolc d2matsolve (T &x, T const &a0, T const &a1, T const &c)

template<typename T >

boolc d2matsolve3D (T &x, T const &a0, T const &a1, T const &c)

When T is a 3D point solving in one plane can fail, hence this routine will attempt to solve in the other two.

boolc d3matsolve (double &x0, double &x1, double &x2, doublec a00, doublec a01, doublec a02, doublec a10, doublec a11, doublec a12, doublec a20, doublec a21, doublec a22, doublec c0, doublec c1, doublec c2)

boolc d3matsolve (point3< double > &x, point3< double > const &a0, point3< double > const &a1, point3< double > const &a2, point3< double > const &c)

template<typename T , typename D >

boolc lineintersection (D &t0, D &t1, T const &a, T const &b, T const &c, T const &d)

boolc solvequadratic (double &t0, double &t1, doublec a, doublec b, doublec c)

a*t^2+b*t+c=0, solve for t.

void matrixmult (point3< double > &y, doublec *m, point3< double > const &x)

void matrixmult (point3< double > &y, point3< double > const &r0, point3< double > const &r1, point3< double > const &r2, point3< double > const &x)

boolc lineSegmentIntersection (point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)

No division operation, return true if the two line segments (p1,p2) and (q1,q2) intersect.

boolc lineIntersection (double &tp, double &tq, point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)

Generally returns true unless the lines are parallel.

boolc lineSegmentIntersection (double &tp, double &tq, point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)

Only calculates both tp and tq if the line segments intersect.

Define Documentation

#define PI 3.1415926535897932384626433

Definition at line 18 of file mathlib.h.

Function Documentation

template<typename T >

boolc d2matsolve	(	T &	x,
		T const &	a0,
		T const &	a1,
		T const &	c
	)			`[inline]`

Definition at line 643 of file mathlib.h.

Referenced by triangle< PT, PD >::fermatpoint(), triangle< PT, PD >::gergonnepoint(), triangle< PT, PD >::incenter(), halfspaceD2< PT, PD >::intersection(), halfspaceD2< PT, PD >::intersectionIn_t(), lineintersection(), triangle< PT, PD >::napoleanpoint(), triangle< PT, PD >::orthocenter(), triangle< PT, PD >::orthocenteri(), test07(), and tetrahedron< PT, PD >::triangleverticiespoint().

00649 {
00650   doublec det = a0.x*a1.y - a1.x*a0.y;
00651   if (det==0.0)
00652     return false;
00653 
00654   doublec detInv = 1.0/det;
00655 
00656   x.x = (c.x*a1.y - c.y*a1.x)*detInv;
00657   x.y = (c.y*a0.x - c.x*a0.y)*detInv;
00658 
00659   return true;
00660 }

template<typename T >

boolc d2matsolve3D	(	T &	x,
		T const &	a0,
		T const &	a1,
		T const &	c
	)			`[inline]`

When T is a 3D point solving in one plane can fail, hence this routine will attempt to solve in the other two.

Definition at line 664 of file mathlib.h.

Referenced by triangle3D< PT, PD >::fermatpoint(), triangle3D< PT, PD >::gergonnepoint(), triangle3D< PT, PD >::incenter(), triangle3D< PT, PD >::napoleanpoint(), triangle3D< PT, PD >::orthocenter(), and triangle3D< PT, PD >::orthocenteri().

00670 {
00671   double det = a0.x*a1.y - a1.x*a0.y;
00672   // <TODO> make this test dependent on number.
00673   if (det!=0.0)
00674   {
00675     doublec detInv = 1.0/det;
00676 
00677     x.x = (c.x*a1.y - c.y*a1.x)*detInv;
00678     x.y = (c.y*a0.x - c.x*a0.y)*detInv;
00679 
00680     return true;
00681   }
00682 
00683   det = a0.x*a1.z - a1.x*a0.z;
00684   if (det!=0.0)
00685   {
00686     doublec detInv = 1.0/det;
00687 
00688     x.x = (c.x*a1.z - c.z*a1.x)*detInv;
00689     x.y = (c.z*a0.x - c.x*a0.z)*detInv;
00690 
00691     return true;
00692   }
00693 
00694   det = a0.y*a1.z - a1.y*a0.z;
00695   if (det!=0.0)
00696   {
00697     doublec detInv = 1.0/det;
00698 
00699     x.x = (c.y*a1.z - c.z*a1.y)*detInv;
00700     x.y = (c.z*a0.y - c.x*a0.z)*detInv;
00701 
00702     return true;
00703   }
00704 
00705   return false;
00706 }

boolc d3matsolve	(	point3< double > &	x,
		point3< double > const &	a0,
		point3< double > const &	a1,
		point3< double > const &	a2,
		point3< double > const &	c
	)

Definition at line 126 of file mathlib.cpp.

References d3matsolve(), point3< T >::x, point3< T >::y, and point3< T >::z.

00133 {
00134   return
00135   d3matsolve
00136   (
00137     x.x,x.y,x.z,
00138     a0.x,a1.x,a2.x,
00139     a0.y,a1.y,a2.y,
00140     a0.z,a1.z,a2.z,
00141     c.x,c.y,c.z
00142   );
00143 }

boolc d3matsolve	(	double &	x0,
		double &	x1,
		double &	x2,
		doublec	a00,
		doublec	a01,
		doublec	a02,
		doublec	a10,
		doublec	a11,
		doublec	a12,
		doublec	a20,
		doublec	a21,
		doublec	a22,
		doublec	c0,
		doublec	c1,
		doublec	c2
	)

Definition at line 102 of file mathlib.cpp.

Referenced by d3matsolve().

00119 {
00120   assert(false);
00121   return false;
00122 }

void integersetdiff	(	vector< uint > &	vecnot,
		vector< uint > &	vec,
		uintc	N
	)

Given the integer set 0,1,.

..,N-1 find vecnot such that all the indexes not in vec are in vecnot. Solve for vecnot.

Definition at line 448 of file mathlib.cpp.

Referenced by mathlibtest::test03().

00453 {
00454   vecnot.clear();
00455   sort(vec.begin(),vec.end());
00456 
00457   uint k=0;
00458   for (uint i=0; i<N; ++i)
00459   {
00460     if (vec[k]==i)
00461       ++k;
00462     else
00463       vecnot.push_back(i);
00464   }
00465 }

boolc lineIntersection	(	double &	tp,
		double &	tq,
		point2< double > const &	p1,
		point2< double > const &	p2,
		point2< double > const &	q1,
		point2< double > const &	q2,
		doublec	zero
	)

Generally returns true unless the lines are parallel.

Only one division performed.

Definition at line 291 of file mathlib.cpp.

References point2< T >::dot(), point2< T >::x, and point2< T >::y.

Referenced by test016(), and test018().

00300 {
00301   point2<double> const a(p2-p1);
00302   point2<double> const aInv(-a.y,a.x);
00303   point2<double> const b(q2-q1);
00304   point2<double> const bInv(-b.y,b.x);
00305 
00306   doublec alpha = b.dot(aInv);
00307 
00308   // Zero test on alpha.
00309   if (alpha + zero > 0.0)
00310   {
00311     if (alpha < zero)
00312       return false;
00313   }
00314 
00315   doublec c = 1.0/alpha;
00316 
00317   point2<double> const w(p1-q1);
00318   tq = c*w.dot(aInv);
00319   tp = c*w.dot(bInv);
00320 
00321   return true;
00322 }

template<typename T , typename D >

boolc lineintersection	(	D &	t0,
		D &	t1,
		T const &	a,
		T const &	b,
		T const &	c,
		T const &	d
	)			`[inline]`

Definition at line 712 of file mathlib.h.

References d2matsolve().

Referenced by d3meshpartition::intersection(), test018(), and test12().

00720 {
00721   T A(b-a);
00722   T B(c-d);
00723   T C(c-a);
00724 
00725   T x;
00726   bool res = d2matsolve(x,A,B,C);
00727   t0 = x.x;
00728   t1 = x.y;
00729 
00730   return res;
00731 }

boolc lineSegmentIntersection	(	double &	tp,
		double &	tq,
		point2< double > const &	p1,
		point2< double > const &	p2,
		point2< double > const &	q1,
		point2< double > const &	q2,
		doublec	zero
	)

Only calculates both tp and tq if the line segments intersect.

tp defines the intersection point p1 + (p2-p1)*tp, similarly tq. Only one division performed.

Definition at line 328 of file mathlib.cpp.

References point2< T >::dot(), point2< T >::x, and point2< T >::y.

00337 {
00338   point2<double> const a(p2-p1);
00339   point2<double> const aInv(-a.y,a.x);
00340   point2<double> const b(q2-q1);
00341   point2<double> const bInv(-b.y,b.x);
00342 
00343   doublec alpha = b.dot(aInv);
00344 
00345   // Zero test on alpha.
00346   if (alpha + zero > 0.0)
00347   {
00348     if (alpha < zero)
00349       return false;
00350   }
00351 
00352   doublec c = 1.0/alpha;
00353 
00354   point2<double> const w(p1-q1);
00355   tq = c*w.dot(aInv);
00356 
00357   if (tq<0.0)
00358     return false;
00359 
00360   if (tq>1.0)
00361     return false;
00362 
00363   tp = c*w.dot(bInv);
00364 
00365   if (tp<0.0)
00366     return false;
00367 
00368   if (tp>1.0)
00369     return false;
00370 
00371   return true;
00372 }

boolc lineSegmentIntersection	(	point2< double > const &	p1,
		point2< double > const &	p2,
		point2< double > const &	q1,
		point2< double > const &	q2,
		doublec	zero
	)

No division operation, return true if the two line segments (p1,p2) and (q1,q2) intersect.

Definition at line 225 of file mathlib.cpp.

References point2< T >::dot(), point2< T >::x, and point2< T >::y.

Referenced by d2linesegment::intersects(), test016(), and test017().

00232 {
00233   point2<double> const a(p2-p1);
00234   point2<double> const aInv(-a.y,a.x);
00235 
00236 //  point2<double> const aInv(p1.y-p2.y,p2.x-p1.x);
00237 
00238 //  point2<double> const b(q2-q1);
00239 //  point2<double> const bInv(b.y*-1.0,b.x);
00240   point2<double> const bInv(q1.y-q2.y,q2.x-q1.x);
00241 
00242   double alpha = a.dot(bInv);
00243 
00244   // Zero test on alpha.
00245   if (alpha + zero > 0.0)
00246   {
00247     if (alpha < zero)
00248       return false;
00249   }
00250 
00251   point2<double> const k(p1-q1);
00252 
00253   double w;
00254 
00255   if (alpha<0.0)
00256   {
00257     alpha *= -1.0;
00258 
00259     w = k.dot(bInv);
00260     if (w<0.0)
00261       return false;
00262     if (w>alpha)
00263       return false;
00264 
00265     w = k.dot(aInv);
00266     if (w<0.0)
00267       return false;
00268 
00269     if (w>alpha)
00270       return false;
00271 
00272     return true;
00273   }
00274 
00275   w = k.dot(bInv)*-1.0;
00276   if (w<0.0)
00277     return false;
00278   if (w>alpha)
00279     return false;
00280 
00281   w = k.dot(aInv)*-1.0;
00282   if (w<0.0)
00283     return false;
00284   if (w>alpha)
00285     return false;
00286 
00287   return true; 
00288 }

void mathcombination	(	double &	result,
		uintc	n,
		uintc	k
	)

n!/((n-k)!k!)

Definition at line 428 of file mathlib.cpp.

Referenced by bernsteinPoly::construct(), and test022().

00433 {
00434 // (n,k) = prod(i=0..n-k-1, (n-i)/(n-k-i) )
00435   result = 1.0;
00436 
00437   uintc ndiffk = n-k;
00438   for (uint i=0; i<ndiffk; ++i)
00439   {
00440     result *= (n-i);
00441     result /= (ndiffk-i);
00442   }
00443 }

void matrixmult	(	point3< double > &	y,
		point3< double > const &	r0,
		point3< double > const &	r1,
		point3< double > const &	r2,
		point3< double > const &	x
	)

Definition at line 210 of file mathlib.cpp.

References point3< T >::dot(), point3< T >::x, point3< T >::y, and point3< T >::z.

00217 {
00218   y.x = r0.dot(x);
00219   y.y = r1.dot(x);
00220   y.z = r2.dot(x);
00221 }

void matrixmult	(	point3< double > &	y,
		doublec *	m,
		point3< double > const &	x
	)

Definition at line 198 of file mathlib.cpp.

References point3< T >::x, point3< T >::y, and point3< T >::z.

Referenced by test014().

00203 {
00204   y.x = m[0]*x.x + m[1]*x.y + m[2]*x.z;
00205   y.y = m[3]*x.x + m[4]*x.y + m[5]*x.z;
00206   y.z = m[6]*x.x + m[7]*x.y + m[8]*x.z;
00207 }

void polygonconvexarea	(	double &	area,
		vector< point2< double > > const &	v
	)

Calculate the area of the convex polygon.

The points are ordered consecutively.

Definition at line 409 of file mathlib.cpp.

00413 {
00414   uintc n = v.size();
00415 
00416   area = v[n-1].x*v[0].y - v[0].x*v[n-1].y;
00417   for (uint i=0; i<n-1; ++i)
00418   {
00419     area += v[i].x*v[i+1].y - v[i+1].x*v[i].y;
00420   }
00421   area *= 0.5;
00422 }

boolc solvequadratic	(	double &	t0,
		double &	t1,
		doublec	a,
		doublec	b,
		doublec	c
	)

a*t^2+b*t+c=0, solve for t.

Definition at line 149 of file mathlib.cpp.

Referenced by tetrahedron< PT, PD >::equilaterali(), d3circlepartition::intersection(), circleLine::intersection2D(), and test11().

00156 {
00157   doublec det = b*b - 4.0 * a * c;
00158   if (det<0.0)
00159     return false;
00160 
00161   double y = sqrt(det);
00162   double a2 = 0.5/a;
00163 
00164   t0 = (-b-y)*a2;
00165   t1 = (-b+y)*a2;
00166 
00167   return true;
00168 }

void tetrahedronvolume	(	double &	vol,
		point3< double > const &	p0,
		point3< double > const &	p1,
		point3< double > const &	p2,
		point3< double > const &	p3
	)

Calculate the volume of the tetrahedron.

Definition at line 375 of file mathlib.cpp.

References crossproduct::eval().

Referenced by test021().

00382 {
00383   //point3<double> a(p1-p0);
00384   //point3<double> b(p2-p0);
00385   //point3<double> c(p3-p0);
00386 
00387 /*
00388 cout << SHOW(p0) << endl;
00389 cout << SHOW(p1) << endl;
00390 cout << SHOW(p2) << endl;
00391 cout << SHOW(p3) << endl;
00392 
00393 cout << SHOW(a) << endl;
00394 cout << SHOW(b) << endl;
00395 cout << SHOW(c) << endl;
00396 */
00397 
00398   point3<double> q;
00399   crossproduct::eval(q,p1-p0,p2-p0);
00400   //crossprod< point3<double> >()(q,p1-p0,p2-p0);
00401 //cout << SHOW(q) << endl;
00402   vol = (p3-p0).dot(q)/6.0;
00403 }

void trianglearea	(	double &	f,
		point3< double > const &	p0,
		point3< double > const &	p1,
		point3< double > const &	p2
	)

Calcluate the triangles area given its three points.

Definition at line 32 of file mathlib.cpp.

References point3< T >::crossproduct(), and point3< T >::distance().

00038 {
00039   point3<double> p;
00040   p.crossproduct(p1-p0,p2-p0);
00041   f = p.distance()*0.5;
00042 }

void trianglearea	(	double &	f,
		doublec	x0,
		doublec	y0,
		doublec	x1,
		doublec	y1,
		doublec	x2,
		doublec	y2
	)

Definition at line 13 of file mathlib.cpp.

Referenced by d4minboundary::eval(), test06(), and test10().

00022 {
00023   double b = point2<double>(x0-x1,y0-y1).distance();
00024   double c = point2<double>(x0-x2,y0-y2).distance();
00025   double dot = x1*x2+y1*y2;
00026   double z(b*c);
00027   f = 0.5*sqrt(z*z-dot*dot);
00028 }

template<typename T >

void unitbound ( T & val ) [inline]

Mape the variable to be inside [0.0,1.0].

Definition at line 23 of file mathlib.h.

Referenced by diskinttest::keyboard03(), and diskinttest::update03lineline().

00024 {
00025   if (val<(T)0)
00026   {
00027     val = (T)0;
00028     return;
00029   }
00030 
00031   if (val>(T)1.0)
00032     val = (T)1.0;
00033 }


Classes
class	crossproduct
	In 3D given two vectors that make a plane find a third vector at right angles to this plane. More...
class	intervalintersection
	Interval overlap tests. More...
class	solverInconsistent< T >
	Solve 1D and 2D linear equations with one move variable than the number of equations. More...
class	surfaceplane< T >
class	convexity< T >
class	orthoproj< T, V, IP, Dim >
class	determinant
	Determinant calculations. More...
class	solver< T >
class	transrotate2D
class	circleLine
	Circle and line intersection test. More...
Defines
#define	PI 3.1415926535897932384626433
Functions
template<typename T >
void	unitbound (T &val)
	Mape the variable to be inside [0.0,1.0].
void	integersetdiff (vector< uint > &vecnot, vector< uint > &vec, uintc N)
	Given the integer set 0,1,.
void	mathcombination (double &result, uintc n, uintc k)
	n!/((n-k)!k!)
void	tetrahedronvolume (double &vol, point3< double > const &p0, point3< double > const &p1, point3< double > const &p2, point3< double > const &p3)
	Calculate the volume of the tetrahedron.
void	polygonconvexarea (double &area, vector< point2< double > > const &v)
	Calculate the area of the convex polygon.
void	trianglearea (double &f, doublec x0, doublec y0, doublec x1, doublec y1, doublec x2, doublec y2)
void	trianglearea (double &f, point3< double > const &p0, point3< double > const &p1, point3< double > const &p2)
	Calcluate the triangles area given its three points.
template<typename T >
boolc	d2matsolve (T &x, T const &a0, T const &a1, T const &c)
template<typename T >
boolc	d2matsolve3D (T &x, T const &a0, T const &a1, T const &c)
	When T is a 3D point solving in one plane can fail, hence this routine will attempt to solve in the other two.
boolc	d3matsolve (double &x0, double &x1, double &x2, doublec a00, doublec a01, doublec a02, doublec a10, doublec a11, doublec a12, doublec a20, doublec a21, doublec a22, doublec c0, doublec c1, doublec c2)
boolc	d3matsolve (point3< double > &x, point3< double > const &a0, point3< double > const &a1, point3< double > const &a2, point3< double > const &c)
template<typename T , typename D >
boolc	lineintersection (D &t0, D &t1, T const &a, T const &b, T const &c, T const &d)
boolc	solvequadratic (double &t0, double &t1, doublec a, doublec b, doublec c)
	at^2+bt+c=0, solve for t.
void	matrixmult (point3< double > &y, doublec *m, point3< double > const &x)
void	matrixmult (point3< double > &y, point3< double > const &r0, point3< double > const &r1, point3< double > const &r2, point3< double > const &x)
boolc	lineSegmentIntersection (point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)
	No division operation, return true if the two line segments (p1,p2) and (q1,q2) intersect.
boolc	lineIntersection (double &tp, double &tq, point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)
	Generally returns true unless the lines are parallel.
boolc	lineSegmentIntersection (double &tp, double &tq, point2< double > const &p1, point2< double > const &p2, point2< double > const &q1, point2< double > const &q2, doublec zero)
	Only calculates both tp and tq if the line segments intersect.