Index: cpp/frams/util/3d.cpp
===================================================================
--- cpp/frams/util/3d.cpp	(revision 247)
+++ cpp/frams/util/3d.cpp	(revision 255)
@@ -7,205 +7,227 @@
 #include "3d.h"
 
-Pt3D operator+(const Pt3D &p1,const Pt3D &p2) {return Pt3D(p1.x+p2.x,p1.y+p2.y,p1.z+p2.z);}
-Pt3D operator-(const Pt3D &p1,const Pt3D &p2) {return Pt3D(p1.x-p2.x,p1.y-p2.y,p1.z-p2.z);}
-
-Pt3D Pt3D_0(0,0,0);
-
-bool Pt3D::report_errors=true;
+Pt3D operator+(const Pt3D &p1, const Pt3D &p2) { return Pt3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z); }
+Pt3D operator-(const Pt3D &p1, const Pt3D &p2) { return Pt3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z); }
+
+Pt3D Pt3D_0(0, 0, 0);
+
+bool Pt3D::report_errors = true;
 
 double Pt3D::operator()() const
 {
-double q=x*x+y*y+z*z;
-if (q<0) {if (report_errors) FramMessage("Pt3D","operator()","sqrt domain error",3); return 0;}
-return sqrt(q);
+	double q = x*x + y*y + z*z;
+	if (q < 0) { if (report_errors) FramMessage("Pt3D", "operator()", "sqrt domain error", 3); return 0; }
+	return sqrt(q);
 }
 
 bool Pt3D::normalize()
 {
-double len=length();
-if (fabs(len)<1e-50) {if (report_errors) FramMessage("Pt3D","normalize()","vector too small",1); x=1;y=0;z=0; return false;}
-operator/=(len);
-return true;
+	double len = length();
+	if (fabs(len) < 1e-50) { if (report_errors) FramMessage("Pt3D", "normalize()", "vector too small", 1); x = 1; y = 0; z = 0; return false; }
+	operator/=(len);
+	return true;
 }
 
 double Pt3D::distanceTo(const Pt3D& p) const
 {
-return sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)+(z-p.z)*(z-p.z));
+	return sqrt((x - p.x)*(x - p.x) + (y - p.y)*(y - p.y) + (z - p.z)*(z - p.z));
 }
 
 double Pt3D::manhattanDistanceTo(const Pt3D& p) const
 {
-return fabs(x-p.x)+fabs(y-p.y)+fabs(z-p.z);
-}
-
-Orient Orient_1(Pt3D(1,0,0),Pt3D(0,1,0),Pt3D(0,0,1));
+	return fabs(x - p.x) + fabs(y - p.y) + fabs(z - p.z);
+}
+
+Orient Orient_1(Pt3D(1, 0, 0), Pt3D(0, 1, 0), Pt3D(0, 0, 1));
 
 // prosty obrot
-void rotate2D(double k,double &x,double &y)
-{double s=sin(k),c=cos(k);
-double t=c*x-s*y; y=s*x+c*y; x=t;}
-
-void rotate2D(double s,double c,double &x,double &y)
-{double t=c*x-s*y; y=s*x+c*y; x=t;}
-
-int Pt3D::getAngle(double dx,double dy,double &wyn)
-{
-if ((fabs(dx)+fabs(dy))<0.0001) return 0;
-wyn=atan2(dy,dx);
-return 1;
-}
-
-void Pt3D::getAngles(const Pt3D& X,const Pt3D& dir)
-{
-Pt3D	t1(X), t2(dir);
-if (getAngle(t1.x,t1.y,z))
-	{	// nie pionowo
-	rotate2D(-z,t1.x,t1.y);
-	rotate2D(-z,t2.x,t2.y);
-	getAngle(t1.x, t1.z, y);
-	}
-else
-	{	// pionowo
-	z=0;
-	if (t1.z<0)
-		y=-M_PI_2; // w dol
-	else
-		y=M_PI_2; // w gore
-	}
-rotate2D(-y,t2.x,t2.z);
-if (!getAngle(t2.z,-t2.y,x))
-	x=0;// to zly wynik ale dobrego nie ma :P
+void rotate2D(double k, double &x, double &y)
+{
+	double s = sin(k), c = cos(k);
+	double t = c*x - s*y; y = s*x + c*y; x = t;
+}
+
+void rotate2D(double s, double c, double &x, double &y)
+{
+	double t = c*x - s*y; y = s*x + c*y; x = t;
+}
+
+int Pt3D::getAngle(double dx, double dy, double &wyn)
+{
+	if ((fabs(dx) + fabs(dy)) < 0.0001) return 0;
+	wyn = atan2(dy, dx);
+	return 1;
+}
+
+void Pt3D::getAngles(const Pt3D& X, const Pt3D& dir)
+{
+	Pt3D	t1(X), t2(dir);
+	if (getAngle(t1.x, t1.y, z)) // non-vertical
+	{	
+		rotate2D(-z, t1.x, t1.y);
+		rotate2D(-z, t2.x, t2.y);
+		getAngle(t1.x, t1.z, y);
+	}
+	else // vertical
+	{	
+		z = 0;
+		if (t1.z < 0)
+			y = -M_PI_2; // down
+		else
+			y = M_PI_2; // up
+	}
+	rotate2D(-y, t2.x, t2.z);
+	if (!getAngle(t2.z, -t2.y, x))
+		x = 0; // incorrect result, but there is no correct one
 }
 
 void Pt3D::getMin(const Pt3D& p)
 {
-if (p.x<x) x=p.x;
-if (p.y<y) y=p.y;
-if (p.z<z) z=p.z;
+	if (p.x < x) x = p.x;
+	if (p.y < y) y = p.y;
+	if (p.z < z) z = p.z;
 }
 void Pt3D::getMax(const Pt3D& p)
 {
-if (p.x>x) x=p.x;
-if (p.y>y) y=p.y;
-if (p.z>z) z=p.z;
-}
-
-void Pt3D::vectorProduct(const Pt3D& a,const Pt3D& b)
-{
-x=a.y*b.z-a.z*b.y;
-y=a.z*b.x-a.x*b.z;
-z=a.x*b.y-a.y*b.x;
-}
-
-void Orient::lookAt(const Pt3D& X,const Pt3D& dir)
-{
-x=X; x.normalize();
-y.vectorProduct(dir,x);
-z.vectorProduct(x,y);
-if ((!y.normalize()) || (!z.normalize()))
-	{
-	y.x=dir.y; y.y=dir.z; y.z=dir.x;
-	z.vectorProduct(x,y);
-	z.normalize();
-	}
-}
-
-// odleglosc 2d
-double d2(double x,double y)
-{
-double q=x*x+y*y;
-if (q<0) {if (Pt3D::report_errors) FramMessage("","d2()","sqrt domain error",3); return 0;}
-return sqrt(q);
+	if (p.x > x) x = p.x;
+	if (p.y > y) y = p.y;
+	if (p.z > z) z = p.z;
+}
+
+void Pt3D::vectorProduct(const Pt3D& a, const Pt3D& b)
+{
+	x = a.y*b.z - a.z*b.y;
+	y = a.z*b.x - a.x*b.z;
+	z = a.x*b.y - a.y*b.x;
+}
+
+void Orient::lookAt(const Pt3D& X, const Pt3D& dir)
+{
+	x = X; x.normalize();
+	y.vectorProduct(dir, x);
+	z.vectorProduct(x, y);
+	if ((!y.normalize()) || (!z.normalize()))
+		lookAt(X);// dir was (nearly?) parallel, there is no good solution, use the x-only variant
+}
+
+void Orient::lookAt(const Pt3D& X)
+{
+	x = X; x.normalize();
+	// "invent" y vector, not parallel to x
+	double ax = fabs(x.x), ay = fabs(x.y), az = fabs(x.z);
+	// find the smallest component
+	if ((ax <= ay) && (ax <= az)) // x
+	{
+		y.x = 0; y.y = -x.z; y.z = x.y; // (0,-z,y)
+	} 
+	if ((ay <= ax) && (ay <= az)) // y
+	{
+		y.x = -x.z; y.y = 0; y.z = x.x; // (-z,0,x)
+	} 
+	else // z
+	{
+		y.x = -x.y; y.y = x.x; y.z = 0; // (-y,x,0)
+	} 
+	y.normalize();
+	z.vectorProduct(x, y);
+}
+
+// 2D distance
+double d2(double x, double y)
+{
+	double q = x*x + y*y;
+	if (q < 0) { if (Pt3D::report_errors) FramMessage("", "d2()", "sqrt domain error", 3); return 0; }
+	return sqrt(q);
 }
 
 Orient::Orient(const Matrix44& m)
 {
-x.x=m[0];  x.y=m[1];  x.z=m[2];
-y.x=m[4];  y.y=m[5];  y.z=m[6];
-z.x=m[8];  z.y=m[9];  z.z=m[10];
+	x.x = m[0];  x.y = m[1];  x.z = m[2];
+	y.x = m[4];  y.y = m[5];  y.z = m[6];
+	z.x = m[8];  z.y = m[9];  z.z = m[10];
 }
 
 void Orient::operator=(const Pt3D &rot)
 {
-*this=Orient_1;
-rotate(rot);
+	*this = Orient_1;
+	rotate(rot);
 }
 
 void Orient::rotate(const Pt3D &v)
 {
-double s,c;
-if (fabs(v.x)>0.0001)
-	{
-	s=sin(v.x); c=cos(v.x);
-	rotate2D(s,c,x.y,x.z);
-	rotate2D(s,c,y.y,y.z);
-	rotate2D(s,c,z.y,z.z);
-	}
-if (fabs(v.y)>0.0001)
-	{
-	s=sin(v.y); c=cos(v.y);
-	rotate2D(s,c,x.x,x.z);
-	rotate2D(s,c,y.x,y.z);
-	rotate2D(s,c,z.x,z.z);
-	}
-if (fabs(v.z)>0.0001)
-	{
-	s=sin(v.z); c=cos(v.z);
-	rotate2D(s,c,x.x,x.y);
-	rotate2D(s,c,y.x,y.y);
-	rotate2D(s,c,z.x,z.y);
-	}
-}
-
-void Orient::transform(Pt3D& target,const Pt3D &s) const
-{
-target.x=s.x*x.x+s.y*y.x+s.z*z.x;
-target.y=s.x*x.y+s.y*y.y+s.z*z.y;
-target.z=s.x*x.z+s.y*y.z+s.z*z.z;
-}
-
-void Orient::revTransform(Pt3D& target,const Pt3D &s) const
-{
-target.x=s.x*x.x+s.y*x.y+s.z*x.z;
-target.y=s.x*y.x+s.y*y.y+s.z*y.z;
-target.z=s.x*z.x+s.y*z.y+s.z*z.z;
-}
-
-void Orient::transform(Orient& target,const Orient& src) const
-{
-transform(target.x,src.x);
-transform(target.y,src.y);
-transform(target.z,src.z);
-}
-
-void Orient::revTransform(Orient& target,const Orient& src) const
-{
-revTransform(target.x,src.x);
-revTransform(target.y,src.y);
-revTransform(target.z,src.z);
+	double s, c;
+	if (fabs(v.x) > 0.0001)
+	{
+		s = sin(v.x); c = cos(v.x);
+		rotate2D(s, c, x.y, x.z);
+		rotate2D(s, c, y.y, y.z);
+		rotate2D(s, c, z.y, z.z);
+	}
+	if (fabs(v.y) > 0.0001)
+	{
+		s = sin(v.y); c = cos(v.y);
+		rotate2D(s, c, x.x, x.z);
+		rotate2D(s, c, y.x, y.z);
+		rotate2D(s, c, z.x, z.z);
+	}
+	if (fabs(v.z) > 0.0001)
+	{
+		s = sin(v.z); c = cos(v.z);
+		rotate2D(s, c, x.x, x.y);
+		rotate2D(s, c, y.x, y.y);
+		rotate2D(s, c, z.x, z.y);
+	}
+}
+
+void Orient::transform(Pt3D& target, const Pt3D &s) const
+{
+	target.x = s.x*x.x + s.y*y.x + s.z*z.x;
+	target.y = s.x*x.y + s.y*y.y + s.z*z.y;
+	target.z = s.x*x.z + s.y*y.z + s.z*z.z;
+}
+
+void Orient::revTransform(Pt3D& target, const Pt3D &s) const
+{
+	target.x = s.x*x.x + s.y*x.y + s.z*x.z;
+	target.y = s.x*y.x + s.y*y.y + s.z*y.z;
+	target.z = s.x*z.x + s.y*z.y + s.z*z.z;
+}
+
+void Orient::transform(Orient& target, const Orient& src) const
+{
+	transform(target.x, src.x);
+	transform(target.y, src.y);
+	transform(target.z, src.z);
+}
+
+void Orient::revTransform(Orient& target, const Orient& src) const
+{
+	revTransform(target.x, src.x);
+	revTransform(target.y, src.y);
+	revTransform(target.z, src.z);
 }
 
 void Orient::getAngles(Pt3D &angles) const
 {
-angles.getAngles(x,z);
+	angles.getAngles(x, z);
 }
 
 bool Orient::normalize()
 {
-bool ret=1;
-y.vectorProduct(z,x);
-z.vectorProduct(x,y);
-if (!x.normalize()) ret=0;
-if (!z.normalize()) ret=0;
-if (!y.normalize()) ret=0;
-return ret;
+	bool ret = 1;
+	y.vectorProduct(z, x);
+	z.vectorProduct(x, y);
+	if (!x.normalize()) ret = 0;
+	if (!z.normalize()) ret = 0;
+	if (!y.normalize()) ret = 0;
+	return ret;
 }
 
 Matrix44::Matrix44(const Orient &rot)
 {
-m[0]=rot.x.x;  m[1]=rot.x.y;  m[2]=rot.x.z;  m[3]=0;
-m[4]=rot.y.x;  m[5]=rot.y.y;  m[6]=rot.y.z;  m[7]=0;
-m[8]=rot.z.x;  m[9]=rot.z.y;  m[10]=rot.z.z; m[11]=0;
-m[12]=0;       m[13]=0;       m[14]=0;       m[15]=1;
+	m[0] = rot.x.x;  m[1] = rot.x.y;  m[2] = rot.x.z;  m[3] = 0;
+	m[4] = rot.y.x;  m[5] = rot.y.y;  m[6] = rot.y.z;  m[7] = 0;
+	m[8] = rot.z.x;  m[9] = rot.z.y;  m[10] = rot.z.z; m[11] = 0;
+	m[12] = 0;       m[13] = 0;       m[14] = 0;       m[15] = 1;
 }
 
Index: cpp/frams/util/3d.h
===================================================================
--- cpp/frams/util/3d.h	(revision 247)
+++ cpp/frams/util/3d.h	(revision 255)
@@ -15,93 +15,95 @@
 \file 3d.h 3d.cpp
 
-   basic 3d classes and operators     
+basic 3D classes and operators
 *********************************/
 
-/// point in 3d space
-
+/// point in 3D space
 class Pt3D
 {
-public:	double x,y,z;
-static bool report_errors;
+public:
+	double x, y, z;
+	static bool report_errors;
 
-Pt3D(double _x,double _y,double _z):x(_x),y(_y),z(_z) {} ///< constructor initializing all coords
-Pt3D(double xyz):x(xyz),y(xyz),z(xyz) {} ///< all coords equal 
-Pt3D() {} ///< coords will be not initialized!
-Pt3D(const Pt3D &p):x(p.x),y(p.y),z(p.z) {} ///< copy from another point
-bool 	operator==(const Pt3D& p)	{return (x==p.x)&&(y==p.y)&&(z==p.z);}
-void 	operator+=(const Pt3D& p)	{x+=p.x;y+=p.y;z+=p.z;}
-void 	operator-=(const Pt3D& p)	{x-=p.x;y-=p.y;z-=p.z;}
-void  	operator*=(double d)	{x*=d;y*=d;z*=d;}
-Pt3D  	operator*(const Pt3D &p) const {return Pt3D(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);}
-void 	operator/=(double d)	{x/=d; y/=d; z/=d;}
-//Pt3D 	operator+(const Pt3D& p) const {return Pt3D(x+p.x,y+p.y,z+p.z);}
-//Pt3D 	operator-(const Pt3D& p) const {return Pt3D(x-p.x,y-p.y,z-p.z);}
-Pt3D 	operator-() const {return Pt3D(-x,-y,-z);}
-Pt3D 	operator*(double d) const {return Pt3D(x*d,y*d,z*d);}
-Pt3D 	operator/(double d) const {return Pt3D(x/d,y/d,z/d);}
-bool    allCoordsLowerThan(const Pt3D& p) const {return (x<p.x)&&(y<p.y)&&(z<p.z);}
-bool    allCoordsHigherThan(const Pt3D& p) const {return (x>p.x)&&(y>p.y)&&(z>p.z);}
-void getMin(const Pt3D& p);
-void getMax(const Pt3D& p);
-/** vector length = \f$\sqrt{x^2+y^2+z^2}\f$  */
-double operator()() const;
-/** vector length = \f$\sqrt{x^2+y^2+z^2}\f$  */
-double length() const {return operator()();}
-double length2() const {return x*x+y*y+z*z;}
-double distanceTo(const Pt3D& p) const;
-double manhattanDistanceTo(const Pt3D& p) const;
-/** calculate angle between (0,0)-(dx,dy), @return 1=ok, 0=can't calculate */
-static int getAngle(double dx,double dy,double &angle);
-/** calculate 3 rotation angles translating (1,0,0) into 'X' and (0,0,1) into 'dir' */
-void getAngles(const Pt3D& X,const Pt3D& dir);
-void vectorProduct(const Pt3D& a,const Pt3D& b);
-Pt3D vectorProduct(const Pt3D& p) const {return (*this)*p;}
-Pt3D entrywiseProduct(const Pt3D &p) const {return Pt3D(x*p.x,y*p.y,z*p.z);} ///< also known as Hadamard product or Schur product
-double dotProduct(const Pt3D& p) const {return x*p.x+y*p.y+z*p.z;}
-bool normalize();
+	Pt3D(double _x, double _y, double _z) :x(_x), y(_y), z(_z) {} ///< constructor initializing all coords
+	Pt3D(double xyz) :x(xyz), y(xyz), z(xyz) {} ///< all coords equal 
+	Pt3D() {} ///< coords will be not initialized!
+	Pt3D(const Pt3D &p) :x(p.x), y(p.y), z(p.z) {} ///< copy from another point
+	bool 	operator==(const Pt3D& p)	{ return (x == p.x) && (y == p.y) && (z == p.z); }
+	void 	operator+=(const Pt3D& p)	{ x += p.x; y += p.y; z += p.z; }
+	void 	operator-=(const Pt3D& p)	{ x -= p.x; y -= p.y; z -= p.z; }
+	void  	operator*=(double d)	{ x *= d; y *= d; z *= d; }
+	Pt3D  	operator*(const Pt3D &p) const { return Pt3D(y*p.z - z*p.y, z*p.x - x*p.z, x*p.y - y*p.x); }
+	void 	operator/=(double d)	{ x /= d; y /= d; z /= d; }
+	//Pt3D 	operator+(const Pt3D& p) const {return Pt3D(x+p.x,y+p.y,z+p.z);}
+	//Pt3D 	operator-(const Pt3D& p) const {return Pt3D(x-p.x,y-p.y,z-p.z);}
+	Pt3D 	operator-() const { return Pt3D(-x, -y, -z); }
+	Pt3D 	operator*(double d) const { return Pt3D(x*d, y*d, z*d); }
+	Pt3D 	operator/(double d) const { return Pt3D(x / d, y / d, z / d); }
+	bool    allCoordsLowerThan(const Pt3D& p) const { return (x < p.x) && (y < p.y) && (z<p.z); }
+	bool    allCoordsHigherThan(const Pt3D& p) const { return (x>p.x) && (y > p.y) && (z > p.z); }
+	void getMin(const Pt3D& p);
+	void getMax(const Pt3D& p);
+	/** vector length = \f$\sqrt{x^2+y^2+z^2}\f$  */
+	double operator()() const;
+	/** vector length = \f$\sqrt{x^2+y^2+z^2}\f$  */
+	double length() const { return operator()(); }
+	double length2() const { return x*x + y*y + z*z; }
+	double distanceTo(const Pt3D& p) const;
+	double manhattanDistanceTo(const Pt3D& p) const;
+	/** calculate angle between (0,0)-(dx,dy), @return 1=ok, 0=can't calculate */
+	static int getAngle(double dx, double dy, double &angle);
+	/** calculate 3 rotation angles translating (1,0,0) into 'X' and (0,0,1) into 'dir' */
+	void getAngles(const Pt3D& X, const Pt3D& dir);
+	void vectorProduct(const Pt3D& a, const Pt3D& b);
+	Pt3D vectorProduct(const Pt3D& p) const { return (*this)*p; }
+	Pt3D entrywiseProduct(const Pt3D &p) const { return Pt3D(x*p.x, y*p.y, z*p.z); } ///< also known as Hadamard product or Schur product
+	double dotProduct(const Pt3D& p) const { return x*p.x + y*p.y + z*p.z; }
+	bool normalize();
 };
-Pt3D operator+(const Pt3D &p1,const Pt3D &p2);
-Pt3D operator-(const Pt3D &p1,const Pt3D &p2);
+
+Pt3D operator+(const Pt3D &p1, const Pt3D &p2);
+Pt3D operator-(const Pt3D &p1, const Pt3D &p2);
 
 class Pt3D_DontReportErrors
 {
-bool state;
+	bool state;
 public:
-Pt3D_DontReportErrors() {state=Pt3D::report_errors; Pt3D::report_errors=false;}
-~Pt3D_DontReportErrors() {Pt3D::report_errors=state;}
+	Pt3D_DontReportErrors() { state = Pt3D::report_errors; Pt3D::report_errors = false; }
+	~Pt3D_DontReportErrors() { Pt3D::report_errors = state; }
 };
 
-///  orientation in 3d space = rotation matrix
-
+///  orientation in 3D space = rotation matrix
 class Matrix44;
 
 class Orient
 {
-public: Pt3D x,y,z; ///< 3 vectors (= 3x3 matrix)
+public:
+	Pt3D x, y, z; ///< 3 vectors (= 3x3 matrix)
 
 	Orient() {}
-	Orient(const Orient& src) {x=src.x; y=src.y; z=src.z;}
-	Orient(const Pt3D& a,const Pt3D& b,const Pt3D& c):x(a),y(b),z(c) {}
-//	Orient(const Pt3D& rot) {*this=rot;}
+	Orient(const Orient& src) { x = src.x; y = src.y; z = src.z; }
+	Orient(const Pt3D& a, const Pt3D& b, const Pt3D& c) :x(a), y(b), z(c) {}
+	//	Orient(const Pt3D& rot) {*this=rot;}
 	Orient(const Matrix44& m);
 	void operator=(const Pt3D &rot);
 	void rotate(const Pt3D &); ///< rotate matrix around 3 axes
 
-	void transform(Pt3D &target,const Pt3D &src) const;	///< transform a vector
-	void revTransform(Pt3D &target,const Pt3D &src) const;	///< reverse transform
-	Pt3D transform(const Pt3D &src) const {Pt3D t; transform(t,src); return t;}
-	Pt3D revTransform(const Pt3D &src) const {Pt3D t; revTransform(t,src); return t;}
+	void transform(Pt3D &target, const Pt3D &src) const;	///< transform a vector
+	void revTransform(Pt3D &target, const Pt3D &src) const;	///< reverse transform
+	Pt3D transform(const Pt3D &src) const { Pt3D t; transform(t, src); return t; }
+	Pt3D revTransform(const Pt3D &src) const { Pt3D t; revTransform(t, src); return t; }
 
-	void transform(Orient& target,const Orient& src) const;    ///< transform other orient
-	void revTransform(Orient& target,const Orient& src) const; ///< reverse transform other orient
-	Orient transform(const Orient& src) const {Orient o; transform(o,src); return o;}    ///< transform other orient
-	Orient revTransform(const Orient& src) const {Orient o; revTransform(o,src); return o;} ///< reverse transform other orient
+	void transform(Orient& target, const Orient& src) const;    ///< transform other orient
+	void revTransform(Orient& target, const Orient& src) const; ///< reverse transform other orient
+	Orient transform(const Orient& src) const { Orient o; transform(o, src); return o; }    ///< transform other orient
+	Orient revTransform(const Orient& src) const { Orient o; revTransform(o, src); return o; } ///< reverse transform other orient
 
-	void transformSelf(const Orient &rot) {Orient tmp; rot.transform(tmp,*this); *this=tmp;}
-	void revTransformSelf(const Orient &rot) {Orient tmp; rot.revTransform(tmp,*this); *this=tmp;}
+	void transformSelf(const Orient &rot) { Orient tmp; rot.transform(tmp, *this); *this = tmp; }
+	void revTransformSelf(const Orient &rot) { Orient tmp; rot.revTransform(tmp, *this); *this = tmp; }
 
 	void getAngles(Pt3D &) const; ///< calculate rotation from current matrix
-	Pt3D getAngles() const {Pt3D ret; getAngles(ret); return ret;}; ///< calculate rotation from current matrix
-	void lookAt(const Pt3D &X,const Pt3D &dir); ///< calculate orientation matrix from 2 vectors
+	Pt3D getAngles() const { Pt3D ret; getAngles(ret); return ret; }; ///< calculate rotation from current matrix
+	void lookAt(const Pt3D &X, const Pt3D &dir); ///< calculate orientation matrix from 2 vectors: X becomes (normalized) Orient.x, dir is the preferred "up" direction (Orient.z). Use lookAt(Pt3D) if only X is relevant.
+	void lookAt(const Pt3D &X); ///< calculate orientation matrix from 1 vector, X becomes (normalized) Orient.x, the other coordinates are deterministic but not continuous. Use lookAt(Pt3D,Pt3D) if you need more control.
 
 	bool normalize();
@@ -111,18 +113,18 @@
 {
 public:
-double m[16];
-Matrix44() {}
-Matrix44(const Matrix44& src) {memcpy(m,src.m,sizeof(m));}
-Matrix44(double *srcm) {memcpy(m,srcm,sizeof(m));}
-Matrix44(const Orient &rot);
+	double m[16];
+	Matrix44() {}
+	Matrix44(const Matrix44& src) { memcpy(m, src.m, sizeof(m)); }
+	Matrix44(double *srcm) { memcpy(m, srcm, sizeof(m)); }
+	Matrix44(const Orient &rot);
 
-const double& operator()(int i,int j) const {return m[i+16*j];}
-const double& operator[](int i) const {return m[i];}
-double& operator()(int i,int j) {return m[i+16*j];}
-double& operator[](int i) {return m[i];}
+	const double& operator()(int i, int j) const { return m[i + 16 * j]; }
+	const double& operator[](int i) const { return m[i]; }
+	double& operator()(int i, int j) { return m[i + 16 * j]; }
+	double& operator[](int i) { return m[i]; }
 
-void operator+=(const Pt3D &); ///< translate matrix
-void operator*=(const Pt3D &); ///< scale matrix
-void operator*=(double sc); ///< scale matrix
+	void operator+=(const Pt3D &); ///< translate matrix
+	void operator*=(const Pt3D &); ///< scale matrix
+	void operator*=(double sc); ///< scale matrix
 };
 
@@ -131,7 +133,7 @@
 extern Matrix44 Matrix44_1; ///< standard unit matrix: 1000 0100 0010 0001
 
-void rotate2D(double,double &,double &); ///< rotate 2d vector, given angle
-void rotate2D(double,double,double &,double &); ///< rotate 2d vector, given sin and cos
-double d2(double,double); ///< distance in 2D
+void rotate2D(double, double &, double &); ///< rotate 2d vector, given angle
+void rotate2D(double, double, double &, double &); ///< rotate 2d vector, given sin and cos
+double d2(double, double); ///< distance in 2D
 
 #endif
