diff --git a/devguide.txt b/devguide.txt index 379c6c28c..66d38afbf 100644 --- a/devguide.txt +++ b/devguide.txt @@ -106,18 +106,9 @@ File Overview Astronomical conversions and a Date class for converting from dd/mm/yyyy to Julian date. -* vecmath.h - Templatized classes for Points, Vectors, and Matrices - -* quaternion.h - Templatized quaternion class - * vecgl.h A handful of template functions which make it easy to use types from - vecmath.h with OpenGL - -* aabox.h - Axis-aligned bounding box class. + Eigen with OpenGL * glext.h * glext.cpp @@ -162,7 +153,7 @@ a pointer to a star from the database. Stars ----- Since there are so many stars in a star database, it's important that the -size of the Star class be kept to a minimum. Thus, in order to avoid the +size of the Star class be kept to a minimum. Thus, in order to avoid the 4 byte vtable overhead Star is not a derived class, even though it might be somewhat useful to have a common base class for stars and planets. All stars have a unique catalog number, a position, a stellar class, and a diff --git a/src/celengine/vertexbuf.h b/src/celengine/vertexbuf.h deleted file mode 100644 index 9a46ddb9c..000000000 --- a/src/celengine/vertexbuf.h +++ /dev/null @@ -1,72 +0,0 @@ -// vertexbuf.h -// -// Copyright (C) 2001, Chris Laurel -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -#ifndef _VERTEXBUF_H_ -#define _VERTEXBUF_H_ - -#include -#include -#include - - -class VertexBuffer -{ - public: - enum { - VertexNormal = 0x01, - VertexColor = 0x02, - VertexColor0 = 0x02, - VertexColor1 = 0x04, - TexCoord0 = 0x08, - TexCoord1 = 0x10, - }; - - class Vertex - { - public: - Point3f point; - Vec3f normal; - Color color; - Point2f texCoords[2]; - }; - - union VertexPart - { - float f; - unsigned char c[4]; - }; - - public: - VertexList(uint32_t _parts, uint32_t initialVertexPoolSize = 0); - ~VertexList(); - - void addVertex(const Vertex& v); - - Color getDiffuseColor() const; - void setDiffuseColor(Color); - - void render(); - - AxisAlignedBox getBoundingBox() const; - void transform(Vec3f translation, float scale); - - private: - uint32_t parts; - uint32_t vertexSize; - - uint32_t nVertices; - uint32_t maxVertices; - VertexPart* vertices; - - Color diffuseColor; - - AxisAlignedBox bbox; -}; - -#endif // _VERTEXBUF_H_ diff --git a/src/celephem/precession.cpp b/src/celephem/precession.cpp index edb887909..e30f19860 100644 --- a/src/celephem/precession.cpp +++ b/src/celephem/precession.cpp @@ -283,94 +283,3 @@ astro::PrecObliquity_P03(double T) return prec; } - - -#define TEST 0 - -#if TEST - -#include - - -int main(int argc, char* argv[]) -{ - double step = 10.0; - - for (int i = -100; i <= 100; i++) - { - // Get time in Julian centuries from J2000 - double T = (double) i * step / 100; - - astro::EclipticPole ecl = astro::EclipticPrecession_P03LP(T); - astro::EclipticPole eclP03 = astro::EclipticPrecession_P03(T); - astro::EclipticAngles eclAnglesP03 = astro::EclipticPrecessionAngles_P03(T); - - //clog << ecl.PA - eclP03.PA << ", " << ecl.QA - eclP03.QA << endl; - ecl = eclP03; - Vec3d v(0.0, 0.0, 0.0); - v.x = ecl.PA * 2.0 * PI / 1296000; - v.y = -ecl.QA * 2.0 * PI / 1296000; - v.z = sqrt(1.0 - v.x * v.x - v.y * v.y); - - Quatd eclRotation = Quatd::vecToVecRotation(Vec3d(0.0, 0.0, 1.0), v); - Quatd eclRotation2 = Quatd::zrotation(-degToRad(eclAnglesP03.PiA / 3600.0)) * - Quatd::xrotation(degToRad(eclAnglesP03.piA / 3600.0)) * - Quatd::zrotation(degToRad(eclAnglesP03.PiA / 3600.0)); - Quatd eclRotation3; - { - // Calculate the angles pi and Pi from the ecliptic pole coordinates - // P and Q: - // P = sin(pi)*sin(Pi) - // Q = sin(pi)*cos(Pi) - double P = eclP03.PA * 2.0 * PI / 1296000; - double Q = eclP03.QA * 2.0 * PI / 1296000; - double piA = asin(sqrt(P * P + Q * Q)); - double PiA = atan2(P, Q); - - // Calculate the rotation from the J2000 ecliptic to the ecliptic - // of date. - eclRotation3 = Quatd::zrotation(-PiA) * - Quatd::xrotation(piA) * - Quatd::zrotation(PiA); - - piA = radToDeg(piA) * 3600; - PiA = radToDeg(PiA) * 3600; - clog << "Pi: " << PiA << ", " << eclAnglesP03.PiA << endl; - clog << "pi: " << piA << ", " << eclAnglesP03.piA << endl; - } - - astro::PrecessionAngles prec = astro::PrecObliquity_P03LP(T); - Quatd p03lpRot = Quatd::xrotation(degToRad(prec.epsA / 3600.0)) * - Quatd::zrotation(-degToRad(prec.pA / 3600.0)); - p03lpRot = p03lpRot * ~eclRotation3; - - astro::PrecessionAngles prec2 = astro::PrecObliquity_P03(T); - //clog << prec.epsA - prec2.epsA << ", " << prec.pA - prec2.pA << endl; - - astro::EquatorialPrecessionAngles precP03 = astro::EquatorialPrecessionAngles_P03(T); - Quatd p03Rot = Quatd::zrotation(-degToRad(precP03.zA / 3600.0)) * - Quatd::yrotation( degToRad(precP03.thetaA / 3600.0)) * - Quatd::zrotation(-degToRad(precP03.zetaA / 3600.0)); - p03Rot = p03Rot * Quatd::xrotation(degToRad(eps0 / 3600.0)); - - Vec3d xaxis(0, 0, 1); - Vec3d v0 = xaxis * p03lpRot.toMatrix3(); - Vec3d v1 = xaxis * p03Rot.toMatrix3(); - - // Show the angle between the J2000 ecliptic pole and the ecliptic - // pole at T centuries from J2000 - double theta0 = acos(xaxis * v0); - double theta1 = acos(xaxis * v1); - - clog << T * 100 << ": " - << radToDeg(acos(v0 * v1)) * 3600 << ", " - << radToDeg(theta0) << ", " - << radToDeg(theta1) << ", " - << radToDeg(theta1 - theta0) * 3600 - << endl; - } - - return 0; -} - -#endif // TEST diff --git a/src/celmath/aabox.h b/src/celmath/aabox.h deleted file mode 100644 index d00f3a377..000000000 --- a/src/celmath/aabox.h +++ /dev/null @@ -1,113 +0,0 @@ -// aabox.h -// -// Copyright (C) 2001, Chris Laurel -// -// Axis-aligned bounding box class -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -#ifndef _AABOX_H_ -#define _AABOX_H_ - -#include - -class AxisAlignedBox -{ - public: - inline AxisAlignedBox(); - AxisAlignedBox(Point3f _min, Point3f _max) : - minimum(_min), maximum(_max) {}; - AxisAlignedBox(Point3f center) : - minimum(center), maximum(center) {}; - - inline Point3f getMinimum() const; - inline Point3f getMaximum() const; - inline Point3f getCenter() const; - inline Vec3f getExtents() const; - - inline bool empty() const; - inline bool contains(const Point3f&) const; - - inline void include(const Point3f&); - inline void include(const AxisAlignedBox&); - - private: - Point3f minimum; - Point3f maximum; -}; - - -AxisAlignedBox::AxisAlignedBox() : - minimum(1.0e20f, 1.0e20f, 1.0e20f), - maximum(-1.0e20f, -1.0e20f, -1.0e20f) -{ -} - -Point3f AxisAlignedBox::getMinimum() const -{ - return minimum; -} - -Point3f AxisAlignedBox::getMaximum() const -{ - return maximum; -} - -Point3f AxisAlignedBox::getCenter() const -{ - return Point3f((minimum.x + maximum.x) * 0.5f, - (minimum.y + maximum.y) * 0.5f, - (minimum.z + maximum.z) * 0.5f); -} - -Vec3f AxisAlignedBox::getExtents() const -{ - return maximum - minimum; -} - -bool AxisAlignedBox::empty() const -{ - return maximum.x < minimum.x || maximum.y < minimum.y || maximum.z < minimum.z; -} - -bool AxisAlignedBox::contains(const Point3f& p) const -{ - return (p.x >= minimum.x && p.x <= maximum.x && - p.y >= minimum.y && p.y <= maximum.y && - p.z >= minimum.z && p.z <= maximum.z); -} - - -void AxisAlignedBox::include(const Point3f& p) -{ - if (p.x < minimum.x) minimum.x = p.x; - if (p.x > maximum.x) maximum.x = p.x; - if (p.y < minimum.y) minimum.y = p.y; - if (p.y > maximum.y) maximum.y = p.y; - if (p.z < minimum.z) minimum.z = p.z; - if (p.z > maximum.z) maximum.z = p.z; -} - -void AxisAlignedBox::include(const AxisAlignedBox& b) -{ - if (b.minimum.x < minimum.x) minimum.x = b.minimum.x; - if (b.maximum.x > maximum.x) maximum.x = b.maximum.x; - if (b.minimum.y < minimum.y) minimum.y = b.minimum.y; - if (b.maximum.y > maximum.y) maximum.y = b.maximum.y; - if (b.minimum.z < minimum.z) minimum.z = b.minimum.z; - if (b.maximum.z > maximum.z) maximum.z = b.maximum.z; -} - -#if 0 -AxisAlignedBox union(const AxisAlignedBox& a, const AxisAlignedBox& b) -{ - AxisAlignedBox box(a); - box.union(b); - return box; -} -#endif - -#endif // _AABOX_H_ diff --git a/src/celmath/capsule.h b/src/celmath/capsule.h deleted file mode 100644 index b2477731b..000000000 --- a/src/celmath/capsule.h +++ /dev/null @@ -1,57 +0,0 @@ -// capsule.h -// -// Implementation of a capsule bounding volume. A capsule is a sphere -// swept along a line, parameterized in this class as a starting point, -// an axis, and a radius. -// -// Copyright (C) 2007, Chris Laurel -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -#ifndef _CELMATH_CAPSULE_H_ -#define _CELMATH_CAPSULE_H_ - -#include "vecmath.h" - -template class Capsule -{ - public: - Capsule(); - Capsule(const Vector3& _axis, T _radius); - Capsule(const Point3& _origin, const Vector3& _axis, T _radius); - - public: - Point3 origin; - Vector3 axis; - T radius; -}; - -typedef Capsule Capsulef; -typedef Capsule Capsuled; - - -template Capsule::Capsule() : - origin(0, 0, 0), axis(0, 0, 0), radius(1) -{ -} - - -template Capsule::Capsule(const Vector3& _axis, - T _radius) : - origin(0, 0, 0), axis(_axis), radius(_radius) -{ -} - - -template Capsule::Capsule(const Point3& _origin, - const Vector3& _axis, - T _radius) : - origin(_origin), axis(_axis), radius(_radius) -{ -} - -#endif // _CELMATH_CAPSULE_H_ - diff --git a/src/celmath/plane.h b/src/celmath/plane.h deleted file mode 100644 index 9666e30fa..000000000 --- a/src/celmath/plane.h +++ /dev/null @@ -1,138 +0,0 @@ -// plane.h -// -// Copyright (C) 2000-2008, Chris Laurel -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -#ifndef _CELMATH_PLANE_H_ -#define _CELMATH_PLANE_H_ - -#include -#include - - -template class Plane -{ - public: - inline Plane(); - inline Plane(const Plane&); - inline Plane(const Vector3&, T); - inline Plane(const Vector3&, const Point3&); - - T distanceTo(const Point3&) const; - T distanceToSegment(const Point3&, const Vector3&) const; - - static Point3 intersection(const Plane&, - const Plane&, - const Plane&); - - public: - Vector3 normal; - T d; -}; - - -typedef Plane Planef; -typedef Plane Planed; - - -template Plane::Plane() : normal(0, 0, 1), d(0) -{ -} - -template Plane::Plane(const Plane& p) : - normal(p.normal), d(p.d) -{ -} - -template Plane::Plane(const Vector3& _normal, T _d) : - normal(_normal), d(_d) -{ -} - -template Plane::Plane(const Vector3& _normal, const Point3& _point) : - normal(_normal) -{ - d = _normal.x * _point.x + _normal.y * _point.y + _normal.z * _point.z; -} - -template T Plane::distanceTo(const Point3& p) const -{ - return normal.x * p.x + normal.y * p.y + normal.z * p.z + d; -} - -// Distance between a plane and a segment defined by orig+dir*t, t <= 0 <= 1 -template T Plane::distanceToSegment(const Point3& origin, - const Vector3& direction) const -{ - T u = (direction * normal); - T dist; - - // Avoid divide by zero; near-zero values shouldn't cause problems - if (u == 0) - { - // All points equidistant; we can just compute distance to origin - dist = distanceTo(origin); - } - else - { - T t = -(d + Vector3(origin.x, origin.y, origin.z) * normal) / u; - if (t < 0) - { - dist = distanceTo(origin); - } - else if (t > 1) - { - dist = distanceTo(origin + direction); - } - else - { - // Segment intersects plane - dist = 0.0; - } - } - - return dist; -} - - -template Plane operator*(const Matrix3& m, const Plane& p) -{ - Vector3 v = m * p.normal; - return Plane(v, p.d); -} - -template Plane operator*(const Plane& p, const Matrix3& m) -{ - Vector3 v = p.normal * m; - return Plane(v, p.d); -} - -template Plane operator*(const Matrix4& m, const Plane& p) -{ - Vector4 v = m * Vector4(p.normal.x, p.normal.y, p.normal.z, p.d); - return Plane(Vector3(v.x, v.y, v.z), v.w); -} - -template Plane operator*(const Plane& p, const Matrix4& m) -{ - Vector4 v = Vector4(p.normal.x, p.normal.y, p.normal.z, p.d) * m; - return Plane(Vector3(v.x, v.y, v.z), v.w); -} - -template Point3 Plane::intersection(const Plane& p0, - const Plane& p1, - const Plane& p2) -{ - T d = Matrix3(p0.normal, p1.normal, p2.normal).determinant(); - - Vector3 v = (p0.d * cross(p1.normal, p2.normal) + - p1.d * cross(p2.normal, p0.normal) + - p2.d * cross(p0.normal, p1.normal)) * (1.0f / d); - return Point3(v.x, v.y, v.z); -} - -#endif // _CELMATH_PLANE_H_ diff --git a/src/celmath/quaternion.h b/src/celmath/quaternion.h deleted file mode 100644 index 693ffd60e..000000000 --- a/src/celmath/quaternion.h +++ /dev/null @@ -1,755 +0,0 @@ -// quaternion.h -// -// Copyright (C) 2000-2006, Chris Laurel -// -// Template-ized quaternion math library. -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -#ifndef _QUATERNION_H_ -#define _QUATERNION_H_ - -#include -#include -#include - - -template class Quat -{ -public: - inline Quat(); - inline Quat(const Quat&); - inline Quat(T); - inline Quat(const Vector3&); - inline Quat(T, const Vector3&); - inline Quat(T, T, T, T); - - inline Quat(const Matrix3&); - - inline Quat& operator+=(Quat); - inline Quat& operator-=(Quat); - inline Quat& operator*=(T); - Quat& operator*=(Quat); - - inline Quat operator~() const; // conjugate - inline Quat operator-() const; - inline Quat operator+() const; - - void setAxisAngle(Vector3 axis, T angle); - - void getAxisAngle(Vector3& axis, T& angle) const; - Matrix4 toMatrix4() const; - Matrix3 toMatrix3() const; - - static Quat slerp(const Quat&, const Quat&, T); - static Quat vecToVecRotation(const Vector3& v0, - const Vector3& v1); - static Quat matrixToQuaternion(const Matrix3& m); - - void rotate(Vector3 axis, T angle); - void xrotate(T angle); - void yrotate(T angle); - void zrotate(T angle); - - bool isPure() const; - bool isReal() const; - T normalize(); - - static Quat xrotation(T); - static Quat yrotation(T); - static Quat zrotation(T); - - static Quat lookAt(const Point3& from, const Point3& to, const Vector3& up); - - T w, x, y, z; -}; - - -typedef Quat Quatf; -typedef Quat Quatd; - - -template Quat::Quat() : w(0), x(0), y(0), z(0) -{ -} - -template Quat::Quat(const Quat& q) : - w(q.w), x(q.x), y(q.y), z(q.z) -{ -} - -template Quat::Quat(T re) : - w(re), x(0), y(0), z(0) -{ -} - -// Create a 'pure' quaternion -template Quat::Quat(const Vector3& im) : - w(0), x(im.x), y(im.y), z(im.z) -{ -} - -template Quat::Quat(T re, const Vector3& im) : - w(re), x(im.x), y(im.y), z(im.z) -{ -} - -template Quat::Quat(T _w, T _x, T _y, T _z) : - w(_w), x(_x), y(_y), z(_z) -{ -} - -// Create a quaternion from a rotation matrix -// TODO: purge this from code--it is replaced by the matrixToQuaternion() -// function. -template Quat::Quat(const Matrix3& m) -{ - T trace = m[0][0] + m[1][1] + m[2][2]; - T root; - - if (trace >= (T) -1.0 + 1.0e-4f) - { - root = (T) sqrt(trace + 1); - w = (T) 0.5 * root; - root = (T) 0.5 / root; - x = (m[2][1] - m[1][2]) * root; - y = (m[0][2] - m[2][0]) * root; - z = (m[1][0] - m[0][1]) * root; - } - else - { - // Identify the largest element of the diagonal - int i = 0; - if (m[1][1] > m[i][i]) - i = 1; - if (m[2][2] > m[i][i]) - i = 2; - int j = (i == 2) ? 0 : i + 1; - int k = (j == 2) ? 0 : j + 1; - - root = (T) sqrt(m[i][i] - m[j][j] - m[k][k] + 1); - T* xyz[3] = { &x, &y, &z }; - *xyz[i] = (T) 0.5 * root; - root = (T) 0.5 / root; - w = (m[k][j] - m[j][k]) * root; - *xyz[j] = (m[j][i] + m[i][j]) * root; - *xyz[k] = (m[k][i] + m[i][k]) * root; - } -} - -template Quat& Quat::operator+=(Quat a) -{ - x += a.x; y += a.y; z += a.z; w += a.w; - return *this; -} - -template Quat& Quat::operator-=(Quat a) -{ - x -= a.x; y -= a.y; z -= a.z; w -= a.w; - return *this; -} - -template Quat& Quat::operator*=(Quat q) -{ - *this = Quat(w * q.w - x * q.x - y * q.y - z * q.z, - w * q.x + x * q.w + y * q.z - z * q.y, - w * q.y + y * q.w + z * q.x - x * q.z, - w * q.z + z * q.w + x * q.y - y * q.x); - - return *this; -} - -template Quat& Quat::operator*=(T s) -{ - x *= s; y *= s; z *= s; w *= s; - return *this; -} - -// conjugate operator -template Quat Quat::operator~() const -{ - return Quat(w, -x, -y, -z); -} - -template Quat Quat::operator-() const -{ - return Quat(-w, -x, -y, -z); -} - -template Quat Quat::operator+() const -{ - return *this; -} - - -template Quat operator+(Quat a, Quat b) -{ - return Quat(a.w + b.w, a.x + b.x, a.y + b.y, a.z + b.z); -} - -template Quat operator-(Quat a, Quat b) -{ - return Quat(a.w - b.w, a.x - b.x, a.y - b.y, a.z - b.z); -} - -template Quat operator*(Quat a, Quat b) -{ - return Quat(a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z, - a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y, - a.w * b.y + a.y * b.w + a.z * b.x - a.x * b.z, - a.w * b.z + a.z * b.w + a.x * b.y - a.y * b.x); -} - -template Quat operator*(T s, Quat q) -{ - return Quat(s * q.w, s * q.x, s * q.y, s * q.z); -} - -template Quat operator*(Quat q, T s) -{ - return Quat(s * q.w, s * q.x, s * q.y, s * q.z); -} - -// equivalent to multiplying by the quaternion (0, v) -template Quat operator*(Vector3 v, Quat q) -{ - return Quat(-v.x * q.x - v.y * q.y - v.z * q.z, - v.x * q.w + v.y * q.z - v.z * q.y, - v.y * q.w + v.z * q.x - v.x * q.z, - v.z * q.w + v.x * q.y - v.y * q.x); -} - -template Quat operator/(Quat q, T s) -{ - return q * (1 / s); -} - -template Quat operator/(Quat a, Quat b) -{ - return a * (~b / abs(b)); -} - - -template bool operator==(Quat a, Quat b) -{ - return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w; -} - -template bool operator!=(Quat a, Quat b) -{ - return a.x != b.x || a.y != b.y || a.z != b.z || a.w != b.w; -} - - -// elementary functions -template Quat conjugate(Quat q) -{ - return Quat(q.w, -q.x, -q.y, -q.z); -} - -template T norm(Quat q) -{ - return q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w; -} - -template T abs(Quat q) -{ - return (T) sqrt(norm(q)); -} - -template Quat exp(Quat q) -{ - if (q.isReal()) - { - return Quat((T) exp(q.w)); - } - else - { - T l = (T) sqrt(q.x * q.x + q.y * q.y + q.z * q.z); - T s = (T) sin(l); - T c = (T) cos(l); - T e = (T) exp(q.w); - T t = e * s / l; - return Quat(e * c, t * q.x, t * q.y, t * q.z); - } -} - -template Quat log(Quat q) -{ - if (q.isReal()) - { - if (q.w > 0) - { - return Quat((T) log(q.w)); - } - else if (q.w < 0) - { - // The log of a negative purely real quaternion has - // infinitely many values, all of the form (ln(-w), PI * I), - // where I is any unit vector. We arbitrarily choose an I - // of (1, 0, 0) here and whereever else a similar choice is - // necessary. Geometrically, the set of roots is a sphere - // of radius PI centered at ln(-w) on the real axis. - return Quat((T) log(-q.w), (T) PI, 0, 0); - } - else - { - // error . . . ln(0) not defined - return Quat(0); - } - } - else - { - T l = (T) sqrt(q.x * q.x + q.y * q.y + q.z * q.z); - T r = (T) sqrt(l * l + q.w * q.w); - T theta = (T) atan2(l, q.w); - T t = theta / l; - return Quat((T) log(r), t * q.x, t * q.y, t * q.z); - } -} - - -template Quat pow(Quat q, T s) -{ - return exp(s * log(q)); -} - - -template Quat pow(Quat q, Quat p) -{ - return exp(p * log(q)); -} - - -template Quat sin(Quat q) -{ - if (q.isReal()) - { - return Quat((T) sin(q.w)); - } - else - { - T l = (T) sqrt(q.x * q.x + q.y * q.y + q.z * q.z); - T m = q.w; - T s = (T) sin(m); - T c = (T) cos(m); - T il = 1 / l; - T e0 = (T) exp(-l); - T e1 = (T) exp(l); - - T c0 = (T) -0.5 * e0 * il * c; - T c1 = (T) 0.5 * e1 * il * c; - - return Quat((T) 0.5 * e0 * s, c0 * q.x, c0 * q.y, c0 * q.z) + - Quat((T) 0.5 * e1 * s, c1 * q.x, c1 * q.y, c1 * q.z); - } -} - -template Quat cos(Quat q) -{ - if (q.isReal()) - { - return Quat((T) cos(q.w)); - } - else - { - T l = (T) sqrt(q.x * q.x + q.y * q.y + q.z * q.z); - T m = q.w; - T s = (T) sin(m); - T c = (T) cos(m); - T il = 1 / l; - T e0 = (T) exp(-l); - T e1 = (T) exp(l); - - T c0 = (T) 0.5 * e0 * il * s; - T c1 = (T) -0.5 * e1 * il * s; - - return Quat((T) 0.5 * e0 * c, c0 * q.x, c0 * q.y, c0 * q.z) + - Quat((T) 0.5 * e1 * c, c1 * q.x, c1 * q.y, c1 * q.z); - } -} - -template Quat sqrt(Quat q) -{ - // In general, the square root of a quaternion has two values, one - // of which is the negative of the other. However, any negative purely - // real quaternion has an infinite number of square roots. - // This function returns the positive root for positive reals and - // the root on the positive i axis for negative reals. - if (q.isReal()) - { - if (q.w >= 0) - return Quat((T) sqrt(q.w), 0, 0, 0); - else - return Quat(0, (T) sqrt(-q.w), 0, 0); - } - else - { - T b = (T) sqrt(q.x * q.x + q.y * q.y + q.z * q.z); - T r = (T) sqrt(q.w * q.w + b * b); - if (q.w >= 0) - { - T m = (T) sqrt((T) 0.5 * (r + q.w)); - T l = b / (2 * m); - T t = l / b; - return Quat(m, q.x * t, q.y * t, q.z * t); - } - else - { - T l = (T) sqrt((T) 0.5 * (r - q.w)); - T m = b / (2 * l); - T t = l / b; - return Quat(m, q.x * t, q.y * t, q.z * t); - } - } -} - -template T real(Quat q) -{ - return q.w; -} - -template Vector3 imag(Quat q) -{ - return Vector3(q.x, q.y, q.z); -} - - -// Quat methods - -template bool Quat::isReal() const -{ - return (x == 0 && y == 0 && z == 0); -} - -template bool Quat::isPure() const -{ - return w == 0; -} - -template T Quat::normalize() -{ - T s = (T) sqrt(w * w + x * x + y * y + z * z); - T invs = (T) 1 / (T) s; - x *= invs; - y *= invs; - z *= invs; - w *= invs; - - return s; -} - -// Set to the unit quaternion representing an axis angle rotation. Assume -// that axis is a unit vector -template void Quat::setAxisAngle(Vector3 axis, T angle) -{ - T s, c; - - Math::sincos(angle * (T) 0.5, s, c); - x = s * axis.x; - y = s * axis.y; - z = s * axis.z; - w = c; -} - - -// Assuming that this a unit quaternion, return the in axis/angle form the -// orientation which it represents. -template void Quat::getAxisAngle(Vector3& axis, - T& angle) const -{ - // The quaternion has the form: - // w = cos(angle/2), (x y z) = sin(angle/2)*axis - - T magSquared = x * x + y * y + z * z; - if (magSquared > (T) 1e-10) - { - T s = (T) 1 / (T) sqrt(magSquared); - axis.x = x * s; - axis.y = y * s; - axis.z = z * s; - if (w <= -1 || w >= 1) - angle = 0; - else - angle = (T) acos(w) * 2; - } - else - { - // The angle is zero, so we pick an arbitrary unit axis - axis.x = 1; - axis.y = 0; - axis.z = 0; - angle = 0; - } -} - - -// Convert this (assumed to be normalized) quaternion to a rotation matrix -template Matrix4 Quat::toMatrix4() const -{ - T wx = w * x * 2; - T wy = w * y * 2; - T wz = w * z * 2; - T xx = x * x * 2; - T xy = x * y * 2; - T xz = x * z * 2; - T yy = y * y * 2; - T yz = y * z * 2; - T zz = z * z * 2; - - return Matrix4(Vector4(1 - yy - zz, xy - wz, xz + wy, 0), - Vector4(xy + wz, 1 - xx - zz, yz - wx, 0), - Vector4(xz - wy, yz + wx, 1 - xx - yy, 0), - Vector4(0, 0, 0, 1)); -} - - -// Convert this (assumed to be normalized) quaternion to a rotation matrix -template Matrix3 Quat::toMatrix3() const -{ - T wx = w * x * 2; - T wy = w * y * 2; - T wz = w * z * 2; - T xx = x * x * 2; - T xy = x * y * 2; - T xz = x * z * 2; - T yy = y * y * 2; - T yz = y * z * 2; - T zz = z * z * 2; - - return Matrix3(Vector3(1 - yy - zz, xy - wz, xz + wy), - Vector3(xy + wz, 1 - xx - zz, yz - wx), - Vector3(xz - wy, yz + wx, 1 - xx - yy)); -} - - -template T dot(Quat a, Quat b) -{ - return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; -} - - -/*! Spherical linear interpolation of two unit quaternions. Designed for - * interpolating rotations, so shortest path between rotations will be - * taken. - */ -template Quat Quat::slerp(const Quat& q0, - const Quat& q1, - T t) -{ - const double Nlerp_Threshold = 0.99999; - - T cosAngle = dot(q0, q1); - - // Assuming the quaternions representat rotations, ensure that we interpolate - // through the shortest path by inverting one of the quaternions if the - // angle between them is negative. - Quat qstart; - if (cosAngle < 0) - { - qstart = -q0; - cosAngle = -cosAngle; - } - else - { - qstart = q0; - } - - // Avoid precision troubles when we're near the limit of acos range and - // perform a linear interpolation followed by a normalize when interpolating - // very small angles. - if (cosAngle > (T) Nlerp_Threshold) - { - Quat q = (1 - t) * qstart + t * q1; - q.normalize(); - return q; - } - - // Below code unnecessary since we've already inverted cosAngle if it's negative. - // It will be necessary if we change slerp to not assume that we want the shortest - // path between two rotations. -#if 0 - // Because of potential rounding errors, we must clamp c to the domain of acos. - if (cosAngle < (T) -1.0) - cosAngle = (T) -1.0; -#endif - - T angle = (T) acos(cosAngle); - T interpolatedAngle = t * angle; - - // qstart and q2 will form an orthonormal basis in the plane of interpolation. - Quat q2 = q1 - qstart * cosAngle; - q2.normalize(); - - return qstart * (T) cos(interpolatedAngle) + q2 * (T) sin(interpolatedAngle); -#if 0 - T s = (T) sin(angle); - T is = (T) 1.0 / s; - - return q0 * ((T) sin((1 - t) * angle) * is) + - q1 * ((T) sin(t * angle) * is); -#endif -} - - -/*! Return a unit quaternion that representing a rotation that will - * rotation v0 to v1 about the axis perpendicular to them both. If the - * vectors point in opposite directions, there is no unique axis and - * (arbitrarily) a rotation about the x axis will be chosen. - */ -template Quat Quat::vecToVecRotation(const Vector3& v0, - const Vector3& v1) -{ - // We need sine and cosine of half the angle between v0 and v1, so - // compute the vector halfway between v0 and v1. The cross product of - // half and v1 gives the imaginary part of the quaternion - // (axis * sin(angle/2)), and the dot product of half and v1 gives - // the real part. - Vector3 half = (v0 + v1) * (T) 0.5; - - T hl = half.length(); - if (hl > (T) 0.0) - { - half = half / hl; // normalize h - - // The magnitude of rotAxis is the sine of half the angle between - // v0 and v1. - Vector3 rotAxis = half ^ v1; - T cosAngle = half * v1; - return Quat(cosAngle, rotAxis.x, rotAxis.y, rotAxis.z); - } - else - { - // The vectors point in exactly opposite directions, so there is - // no unique axis of rotation. Rotating v0 180 degrees about any - // axis will map it to v1; we'll choose the x-axis. - return Quat((T) 0.0, (T) 1.0, (T) 0.0, (T) 0.0); - } -} - - -/*! Create a quaternion from a rotation matrix - */ -template Quat Quat::matrixToQuaternion(const Matrix3& m) -{ - Quat q; - T trace = m[0][0] + m[1][1] + m[2][2]; - T root; - T epsilon = std::numeric_limits::epsilon() * (T) 1e3; - - if (trace >= epsilon - 1) - { - root = (T) sqrt(trace + 1); - q.w = (T) 0.5 * root; - root = (T) 0.5 / root; - q.x = (m[2][1] - m[1][2]) * root; - q.y = (m[0][2] - m[2][0]) * root; - q.z = (m[1][0] - m[0][1]) * root; - } - else - { - // Set i to the largest element of the diagonal - int i = 0; - if (m[1][1] > m[i][i]) - i = 1; - if (m[2][2] > m[i][i]) - i = 2; - int j = (i == 2) ? 0 : i + 1; - int k = (j == 2) ? 0 : j + 1; - - root = (T) sqrt(m[i][i] - m[j][j] - m[k][k] + 1); - T* xyz[3] = { &q.x, &q.y, &q.z }; - *xyz[i] = (T) 0.5 * root; - root = (T) 0.5 / root; - q.w = (m[k][j] - m[j][k]) * root; - *xyz[j] = (m[j][i] + m[i][j]) * root; - *xyz[k] = (m[k][i] + m[i][k]) * root; - } - - return q; -} - - -/*! Assuming that this is a unit quaternion representing an orientation, - * apply a rotation of angle radians about the specfied axis - */ -template void Quat::rotate(Vector3 axis, T angle) -{ - Quat q; - q.setAxisAngle(axis, angle); - *this = q * *this; -} - - -// Assuming that this is a unit quaternion representing an orientation, -// apply a rotation of angle radians about the x-axis -template void Quat::xrotate(T angle) -{ - T s, c; - - Math::sincos(angle * (T) 0.5, s, c); - *this = Quat(c, s, 0, 0) * *this; -} - -// Assuming that this is a unit quaternion representing an orientation, -// apply a rotation of angle radians about the y-axis -template void Quat::yrotate(T angle) -{ - T s, c; - - Math::sincos(angle * (T) 0.5, s, c); - *this = Quat(c, 0, s, 0) * *this; -} - -// Assuming that this is a unit quaternion representing an orientation, -// apply a rotation of angle radians about the z-axis -template void Quat::zrotate(T angle) -{ - T s, c; - - Math::sincos(angle * (T) 0.5, s, c); - *this = Quat(c, 0, 0, s) * *this; -} - - -template Quat Quat::xrotation(T angle) -{ - T s, c; - Math::sincos(angle * (T) 0.5, s, c); - return Quat(c, s, 0, 0); -} - -template Quat Quat::yrotation(T angle) -{ - T s, c; - Math::sincos(angle * (T) 0.5, s, c); - return Quat(c, 0, s, 0); -} - -template Quat Quat::zrotation(T angle) -{ - T s, c; - Math::sincos(angle * (T) 0.5, s, c); - return Quat(c, 0, 0, s); -} - -/*! Determine an orientation that will make the negative z-axis point from - * from the observer to the target, with the y-axis pointing in direction - * of the component of 'up' that is orthogonal to the z-axis. - */ -template Quat -Quat::lookAt(const Point3& from, const Point3& to, const Vector3& up) -{ - Vector3 n = to - from; - n.normalize(); - Vector3 v = n ^ up; - v.normalize(); - Vector3 u = v ^ n; - - return Quat::matrixToQuaternion(Matrix3(v, u, -n)); -} - -#endif // _QUATERNION_H_ diff --git a/src/celmath/vecmath.h b/src/celmath/vecmath.h deleted file mode 100644 index b2fc5ae2d..000000000 --- a/src/celmath/vecmath.h +++ /dev/null @@ -1,1075 +0,0 @@ -// vecmath.h -// -// Copyright (C) 2000, Chris Laurel -// -// Basic templatized vector and matrix math library. -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - - -#ifndef _VECMATH_H_ -#define _VECMATH_H_ - -#include - - -template class Point3; - -template class Vector3 -{ -public: - inline Vector3(); - inline Vector3(const Vector3&); - inline Vector3(T, T, T); - inline Vector3(T*); - - inline T& operator[](int) const; - inline Vector3& operator+=(const Vector3&); - inline Vector3& operator-=(const Vector3&); - inline Vector3& operator*=(T); - - inline Vector3 operator-() const; - inline Vector3 operator+() const; - - inline T length() const; - inline T lengthSquared() const; - inline void normalize(); - - T x, y, z; -}; - -template class Point3 -{ -public: - inline Point3(); - inline Point3(const Point3&); - inline Point3(T, T, T); - inline Point3(T*); - - inline T& operator[](int) const; - inline Point3& operator+=(const Vector3&); - inline Point3& operator-=(const Vector3&); - inline Point3& operator*=(T); - - inline T distanceTo(const Point3&) const; - inline T distanceToSquared(const Point3&) const; - inline T distanceFromOrigin() const; - inline T distanceFromOriginSquared() const; - - T x, y, z; -}; - - -template class Point2; - -template class Vector2 -{ -public: - inline Vector2(); - inline Vector2(T, T); - - T x, y; -}; - -template class Point2 -{ -public: - inline Point2(); - inline Point2(T, T); - - T x, y; -}; - - -template class Vector4 -{ - public: - inline Vector4(); - inline Vector4(T, T, T, T); - - inline T& operator[](int) const; - inline Vector4& operator+=(const Vector4&); - inline Vector4& operator-=(const Vector4&); - inline Vector4& operator*=(T); - - inline Vector4 operator-() const; - inline Vector4 operator+() const; - - T x, y, z, w; -}; - - -template class Matrix4 -{ - public: - Matrix4(); - Matrix4(const Vector4&, const Vector4&, - const Vector4&, const Vector4&); - Matrix4(const Matrix4& m); - - inline const Vector4& operator[](int) const; - inline Vector4 row(int) const; - inline Vector4 column(int) const; - - static Matrix4 identity(); - static Matrix4 translation(const Point3&); - static Matrix4 translation(const Vector3&); - static Matrix4 rotation(const Vector3&, T); - static Matrix4 xrotation(T); - static Matrix4 yrotation(T); - static Matrix4 zrotation(T); - static Matrix4 scaling(const Vector3&); - static Matrix4 scaling(T); - - void translate(const Point3&); - - Matrix4 transpose() const; - Matrix4 inverse() const; - - Vector4 r[4]; -}; - - -template class Matrix3 -{ - public: - Matrix3(); - Matrix3(const Vector3&, const Vector3&, const Vector3&); - template Matrix3(const Matrix3&); - - static Matrix3 xrotation(T); - static Matrix3 yrotation(T); - static Matrix3 zrotation(T); - static Matrix3 scaling(const Vector3&); - static Matrix3 scaling(T); - - inline const Vector3& operator[](int) const; - inline Vector3 row(int) const; - inline Vector3 column(int) const; - - inline Matrix3& operator*=(T); - - static Matrix3 identity(); - - Matrix3 transpose() const; - Matrix3 inverse() const; - T determinant() const; - - // template operator Matrix3() const; - - Vector3 r[3]; -}; - - -typedef Vector3 Vec3f; -typedef Vector3 Vec3d; -typedef Point3 Point3f; -typedef Point3 Point3d; -typedef Vector2 Vec2f; -typedef Point2 Point2f; -typedef Vector4 Vec4f; -typedef Vector4 Vec4d; -typedef Matrix4 Mat4f; -typedef Matrix4 Mat4d; -typedef Matrix3 Mat3f; -typedef Matrix3 Mat3d; - - -//**** Vector3 constructors - -template Vector3::Vector3() : x(0), y(0), z(0) -{ -} - -template Vector3::Vector3(const Vector3& v) : - x(v.x), y(v.y), z(v.z) -{ -} - -template Vector3::Vector3(T _x, T _y, T _z) : x(_x), y(_y), z(_z) -{ -} - -template Vector3::Vector3(T* v) : x(v[0]), y(v[1]), z(v[2]) -{ -} - -//**** Vector3 operators - -template T& Vector3::operator[](int n) const -{ - // Not portable--I'll write a new version when I try to compile on a - // platform where it bombs. - return ((T*) this)[n]; -} - -template Vector3& Vector3::operator+=(const Vector3& a) -{ - x += a.x; y += a.y; z += a.z; - return *this; -} - -template Vector3& Vector3::operator-=(const Vector3& a) -{ - x -= a.x; y -= a.y; z -= a.z; - return *this; -} - -template Vector3& Vector3::operator*=(T s) -{ - x *= s; y *= s; z *= s; - return *this; -} - -template Vector3 Vector3::operator-() const -{ - return Vector3(-x, -y, -z); -} - -template Vector3 Vector3::operator+() const -{ - return *this; -} - - -template Vector3 operator+(const Vector3& a, const Vector3& b) -{ - return Vector3(a.x + b.x, a.y + b.y, a.z + b.z); -} - -template Vector3 operator-(const Vector3& a, const Vector3& b) -{ - return Vector3(a.x - b.x, a.y - b.y, a.z - b.z); -} - -template Vector3 operator*(T s, const Vector3& v) -{ - return Vector3(s * v.x, s * v.y, s * v.z); -} - -template Vector3 operator*(const Vector3& v, T s) -{ - return Vector3(s * v.x, s * v.y, s * v.z); -} - -// dot product -template T operator*(const Vector3& a, const Vector3& b) -{ - return a.x * b.x + a.y * b.y + a.z * b.z; -} - -// cross product -template Vector3 operator^(const Vector3& a, const Vector3& b) -{ - return Vector3(a.y * b.z - a.z * b.y, - a.z * b.x - a.x * b.z, - a.x * b.y - a.y * b.x); -} - -template bool operator==(const Vector3& a, const Vector3& b) -{ - return a.x == b.x && a.y == b.y && a.z == b.z; -} - -template bool operator!=(const Vector3& a, const Vector3& b) -{ - return a.x != b.x || a.y != b.y || a.z != b.z; -} - -template Vector3 operator/(const Vector3& v, T s) -{ - T is = 1 / s; - return Vector3(is * v.x, is * v.y, is * v.z); -} - -template T dot(const Vector3& a, const Vector3& b) -{ - return a.x * b.x + a.y * b.y + a.z * b.z; -} - -template Vector3 cross(const Vector3& a, const Vector3& b) -{ - return Vector3(a.y * b.z - a.z * b.y, - a.z * b.x - a.x * b.z, - a.x * b.y - a.y * b.x); -} - -template T Vector3::length() const -{ - return (T) sqrt(x * x + y * y + z * z); -} - -template T Vector3::lengthSquared() const -{ - return x * x + y * y + z * z; -} - -template void Vector3::normalize() -{ - T s = 1 / (T) sqrt(x * x + y * y + z * z); - x *= s; - y *= s; - z *= s; -} - - -//**** Point3 constructors - -template Point3::Point3() : x(0), y(0), z(0) -{ -} - -template Point3::Point3(const Point3& p) : - x(p.x), y(p.y), z(p.z) -{ -} - -template Point3::Point3(T _x, T _y, T _z) : x(_x), y(_y), z(_z) -{ -} - -template Point3::Point3(T* p) : x(p[0]), y(p[1]), z(p[2]) -{ -} - - -//**** Point3 operators - -template T& Point3::operator[](int n) const -{ - // Not portable--I'll write a new version when I try to compile on a - // platform where it bombs. - return ((T*) this)[n]; -} - -template Vector3 operator-(const Point3& a, const Point3& b) -{ - return Vector3(a.x - b.x, a.y - b.y, a.z - b.z); -} - -template Point3& Point3::operator+=(const Vector3& v) -{ - x += v.x; y += v.y; z += v.z; - return *this; -} - -template Point3& Point3::operator-=(const Vector3& v) -{ - x -= v.x; y -= v.y; z -= v.z; - return *this; -} - -template Point3& Point3::operator*=(T s) -{ - x *= s; y *= s; z *= s; - return *this; -} - -template bool operator==(const Point3& a, const Point3& b) -{ - return a.x == b.x && a.y == b.y && a.z == b.z; -} - -template bool operator!=(const Point3& a, const Point3& b) -{ - return a.x != b.x || a.y != b.y || a.z != b.z; -} - -template Point3 operator+(const Point3& p, const Vector3& v) -{ - return Point3(p.x + v.x, p.y + v.y, p.z + v.z); -} - -template Point3 operator-(const Point3& p, const Vector3& v) -{ - return Point3(p.x - v.x, p.y - v.y, p.z - v.z); -} - -// Naughty naughty . . . remove these since they aren't proper -// point methods--changing the implicit homogenous coordinate isn't -// allowed. -template Point3 operator*(const Point3& p, T s) -{ - return Point3(p.x * s, p.y * s, p.z * s); -} - -template Point3 operator*(T s, const Point3& p) -{ - return Point3(p.x * s, p.y * s, p.z * s); -} - - -//**** Point3 methods - -template T Point3::distanceTo(const Point3& p) const -{ - return (T) sqrt((p.x - x) * (p.x - x) + - (p.y - y) * (p.y - y) + - (p.z - z) * (p.z - z)); -} - -template T Point3::distanceToSquared(const Point3& p) const -{ - return ((p.x - x) * (p.x - x) + - (p.y - y) * (p.y - y) + - (p.z - z) * (p.z - z)); -} - -template T Point3::distanceFromOrigin() const -{ - return (T) sqrt(x * x + y * y + z * z); -} - -template T Point3::distanceFromOriginSquared() const -{ - return x * x + y * y + z * z; -} - - - -//**** Vector2 constructors -template Vector2::Vector2() : x(0), y(0) -{ -} - -template Vector2::Vector2(T _x, T _y) : x(_x), y(_y) -{ -} - - -//**** Vector2 operators -template bool operator==(const Vector2& a, const Vector2& b) -{ - return a.x == b.x && a.y == b.y; -} - -template bool operator!=(const Vector2& a, const Vector2& b) -{ - return a.x != b.x || a.y != b.y; -} - - -//**** Point2 constructors - -template Point2::Point2() : x(0), y(0) -{ -} - -template Point2::Point2(T _x, T _y) : x(_x), y(_y) -{ -} - - -//**** Point2 operators - -template bool operator==(const Point2& a, const Point2& b) -{ - return a.x == b.x && a.y == b.y; -} - -template bool operator!=(const Point2& a, const Point2& b) -{ - return a.x != b.x || a.y != b.y; -} - - -//**** Vector4 constructors - -template Vector4::Vector4() : x(0), y(0), z(0), w(0) -{ -} - -template Vector4::Vector4(T _x, T _y, T _z, T _w) : - x(_x), y(_y), z(_z), w(_w) -{ -} - - -//**** Vector4 operators - -template T& Vector4::operator[](int n) const -{ - // Not portable--I'll write a new version when I try to compile on a - // platform where it bombs. - return ((T*) this)[n]; -} - -template bool operator==(const Vector4& a, const Vector4& b) -{ - return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w; -} - -template bool operator!=(const Vector4& a, const Vector4& b) -{ - return a.x != b.x || a.y != b.y || a.z != b.z || a.w != b.w; -} - -template Vector4& Vector4::operator+=(const Vector4& a) -{ - x += a.x; y += a.y; z += a.z; w += a.w; - return *this; -} - -template Vector4& Vector4::operator-=(const Vector4& a) -{ - x -= a.x; y -= a.y; z -= a.z; w -= a.w; - return *this; -} - -template Vector4& Vector4::operator*=(T s) -{ - x *= s; y *= s; z *= s; w *= s; - return *this; -} - -template Vector4 Vector4::operator-() const -{ - return Vector4(-x, -y, -z, -w); -} - -template Vector4 Vector4::operator+() const -{ - return *this; -} - -template Vector4 operator+(const Vector4& a, const Vector4& b) -{ - return Vector4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); -} - -template Vector4 operator-(const Vector4& a, const Vector4& b) -{ - return Vector4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); -} - -template Vector4 operator*(T s, const Vector4& v) -{ - return Vector4(s * v.x, s * v.y, s * v.z, s * v.w); -} - -template Vector4 operator*(const Vector4& v, T s) -{ - return Vector4(s * v.x, s * v.y, s * v.z, s * v.w); -} - -// dot product -template T operator*(const Vector4& a, const Vector4& b) -{ - return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; -} - -template T dot(const Vector4& a, const Vector4& b) -{ - return a * b; -} - - - -//**** Matrix3 constructors - -template Matrix3::Matrix3() -{ - r[0] = Vector3(0, 0, 0); - r[1] = Vector3(0, 0, 0); - r[2] = Vector3(0, 0, 0); -} - - -template Matrix3::Matrix3(const Vector3& r0, - const Vector3& r1, - const Vector3& r2) -{ - r[0] = r0; - r[1] = r1; - r[2] = r2; -} - - -#if 0 -template Matrix3::Matrix3(const Matrix3& m) -{ -#if 0 - r[0] = m.r[0]; - r[1] = m.r[1]; - r[2] = m.r[2]; -#endif - r[0].x = m.r[0].x; r[0].y = m.r[0].y; r[0].z = m.r[0].z; - r[1].x = m.r[1].x; r[1].y = m.r[1].y; r[1].z = m.r[1].z; - r[2].x = m.r[2].x; r[2].y = m.r[2].y; r[2].z = m.r[2].z; -} -#endif - - -//**** Matrix3 operators - -template const Vector3& Matrix3::operator[](int n) const -{ - // return const_cast&>(r[n]); - return r[n]; -} - -template Vector3 Matrix3::row(int n) const -{ - return r[n]; -} - -template Vector3 Matrix3::column(int n) const -{ - return Vector3(r[0][n], r[1][n], r[2][n]); -} - -template Matrix3& Matrix3::operator*=(T s) -{ - r[0] *= s; - r[1] *= s; - r[2] *= s; - return *this; -} - - -// pre-multiply column vector by a 3x3 matrix -template Vector3 operator*(const Matrix3& m, const Vector3& v) -{ - return Vector3(m.r[0].x * v.x + m.r[0].y * v.y + m.r[0].z * v.z, - m.r[1].x * v.x + m.r[1].y * v.y + m.r[1].z * v.z, - m.r[2].x * v.x + m.r[2].y * v.y + m.r[2].z * v.z); -} - - -// post-multiply row vector by a 3x3 matrix -template Vector3 operator*(const Vector3& v, const Matrix3& m) -{ - return Vector3(m.r[0].x * v.x + m.r[1].x * v.y + m.r[2].x * v.z, - m.r[0].y * v.x + m.r[1].y * v.y + m.r[2].y * v.z, - m.r[0].z * v.x + m.r[1].z * v.y + m.r[2].z * v.z); -} - - -// pre-multiply column point by a 3x3 matrix -template Point3 operator*(const Matrix3& m, const Point3& p) -{ - return Point3(m.r[0].x * p.x + m.r[0].y * p.y + m.r[0].z * p.z, - m.r[1].x * p.x + m.r[1].y * p.y + m.r[1].z * p.z, - m.r[2].x * p.x + m.r[2].y * p.y + m.r[2].z * p.z); -} - - -// post-multiply row point by a 3x3 matrix -template Point3 operator*(const Point3& p, const Matrix3& m) -{ - return Point3(m.r[0].x * p.x + m.r[1].x * p.y + m.r[2].x * p.z, - m.r[0].y * p.x + m.r[1].y * p.y + m.r[2].y * p.z, - m.r[0].z * p.x + m.r[1].z * p.y + m.r[2].z * p.z); -} - - -template Matrix3 operator*(const Matrix3& a, - const Matrix3& b) -{ -#define MATMUL(R, C) (a[R].x * b[0].C + a[R].y * b[1].C + a[R].z * b[2].C) - return Matrix3(Vector3(MATMUL(0, x), MATMUL(0, y), MATMUL(0, z)), - Vector3(MATMUL(1, x), MATMUL(1, y), MATMUL(1, z)), - Vector3(MATMUL(2, x), MATMUL(2, y), MATMUL(2, z))); -#undef MATMUL -} - - -template Matrix3 operator+(const Matrix3& a, - const Matrix3& b) -{ - return Matrix3(a.r[0] + b.r[0], - a.r[1] + b.r[1], - a.r[2] + b.r[2]); -} - - -template Matrix3 Matrix3::identity() -{ - return Matrix3(Vector3(1, 0, 0), - Vector3(0, 1, 0), - Vector3(0, 0, 1)); -} - - -template Matrix3 Matrix3::transpose() const -{ - return Matrix3(Vector3(r[0].x, r[1].x, r[2].x), - Vector3(r[0].y, r[1].y, r[2].y), - Vector3(r[0].z, r[1].z, r[2].z)); -} - - -template T det2x2(T a, T b, T c, T d) -{ - return a * d - b * c; -} - -template T Matrix3::determinant() const -{ - return (r[0].x * r[1].y * r[2].z + - r[0].y * r[1].z * r[2].x + - r[0].z * r[1].x * r[2].y - - r[0].z * r[1].y * r[2].x - - r[0].x * r[1].z * r[2].y - - r[0].y * r[1].x * r[2].z); -} - - -template Matrix3 Matrix3::inverse() const -{ - Matrix3 adjoint; - - // Just use Cramer's rule for now . . . - adjoint.r[0].x = det2x2(r[1].y, r[1].z, r[2].y, r[2].z); - adjoint.r[0].y = -det2x2(r[1].x, r[1].z, r[2].x, r[2].z); - adjoint.r[0].z = det2x2(r[1].x, r[1].y, r[2].x, r[2].y); - adjoint.r[1].x = -det2x2(r[0].y, r[0].z, r[2].y, r[2].z); - adjoint.r[1].y = det2x2(r[0].x, r[0].z, r[2].x, r[2].z); - adjoint.r[1].z = -det2x2(r[0].x, r[0].y, r[2].x, r[2].y); - adjoint.r[2].x = det2x2(r[0].y, r[0].z, r[1].y, r[1].z); - adjoint.r[2].y = -det2x2(r[0].x, r[0].z, r[1].x, r[1].z); - adjoint.r[2].z = det2x2(r[0].x, r[0].y, r[1].x, r[1].y); - adjoint *= 1 / determinant(); - - return adjoint; -} - - -template Matrix3 Matrix3::xrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix3(Vector3(1, 0, 0), - Vector3(0, c, -s), - Vector3(0, s, c)); -} - - -template Matrix3 Matrix3::yrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix3(Vector3(c, 0, s), - Vector3(0, 1, 0), - Vector3(-s, 0, c)); -} - - -template Matrix3 Matrix3::zrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix3(Vector3(c, -s, 0), - Vector3(s, c, 0), - Vector3(0, 0, 1)); -} - - -template Matrix3 Matrix3::scaling(const Vector3& scale) -{ - return Matrix3(Vector3(scale.x, 0, 0), - Vector3(0, scale.y, 0), - Vector3(0, 0, scale.z)); -} - - -template Matrix3 Matrix3::scaling(T scale) -{ - return scaling(Vector3(scale, scale, scale)); -} - - -/*********************************************** - ** - ** Matrix4 methods - ** - ***********************************************/ - -template Matrix4::Matrix4() -{ - r[0] = Vector4(0, 0, 0, 0); - r[1] = Vector4(0, 0, 0, 0); - r[2] = Vector4(0, 0, 0, 0); - r[3] = Vector4(0, 0, 0, 0); -} - - -template Matrix4::Matrix4(const Vector4& v0, - const Vector4& v1, - const Vector4& v2, - const Vector4& v3) -{ - r[0] = v0; - r[1] = v1; - r[2] = v2; - r[3] = v3; -} - - -template Matrix4::Matrix4(const Matrix4& m) -{ - r[0] = m.r[0]; - r[1] = m.r[1]; - r[2] = m.r[2]; - r[3] = m.r[3]; -} - - -template const Vector4& Matrix4::operator[](int n) const -{ - return r[n]; - // return const_cast&>(r[n]); -} - -template Vector4 Matrix4::row(int n) const -{ - return r[n]; -} - -template Vector4 Matrix4::column(int n) const -{ - return Vector4(r[0][n], r[1][n], r[2][n], r[3][n]); -} - - -template Matrix4 Matrix4::identity() -{ - return Matrix4(Vector4(1, 0, 0, 0), - Vector4(0, 1, 0, 0), - Vector4(0, 0, 1, 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::translation(const Point3& p) -{ - return Matrix4(Vector4(1, 0, 0, 0), - Vector4(0, 1, 0, 0), - Vector4(0, 0, 1, 0), - Vector4(p.x, p.y, p.z, 1)); -} - - -template Matrix4 Matrix4::translation(const Vector3& v) -{ - return Matrix4(Vector4(1, 0, 0, 0), - Vector4(0, 1, 0, 0), - Vector4(0, 0, 1, 0), - Vector4(v.x, v.y, v.z, 1)); -} - - -template void Matrix4::translate(const Point3& p) -{ - r[3].x += p.x; - r[3].y += p.y; - r[3].z += p.z; -} - - -template Matrix4 Matrix4::rotation(const Vector3& axis, - T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - T t = 1 - c; - - return Matrix4(Vector4(t * axis.x * axis.x + c, - t * axis.x * axis.y - s * axis.z, - t * axis.x * axis.z + s * axis.y, - 0), - Vector4(t * axis.x * axis.y + s * axis.z, - t * axis.y * axis.y + c, - t * axis.y * axis.z - s * axis.x, - 0), - Vector4(t * axis.x * axis.z - s * axis.y, - t * axis.y * axis.z + s * axis.x, - t * axis.z * axis.z + c, - 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::xrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix4(Vector4(1, 0, 0, 0), - Vector4(0, c, -s, 0), - Vector4(0, s, c, 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::yrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix4(Vector4(c, 0, s, 0), - Vector4(0, 1, 0, 0), - Vector4(-s, 0, c, 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::zrotation(T angle) -{ - T c = (T) cos(angle); - T s = (T) sin(angle); - - return Matrix4(Vector4(c, -s, 0, 0), - Vector4(s, c, 0, 0), - Vector4(0, 0, 1, 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::scaling(const Vector3& scale) -{ - return Matrix4(Vector4(scale.x, 0, 0, 0), - Vector4(0, scale.y, 0, 0), - Vector4(0, 0, scale.z, 0), - Vector4(0, 0, 0, 1)); -} - - -template Matrix4 Matrix4::scaling(T scale) -{ - return scaling(Vector3(scale, scale, scale)); -} - - -// multiply column vector by a 4x4 matrix -template Vector3 operator*(const Matrix4& m, const Vector3& v) -{ - return Vector3(m.r[0].x * v.x + m.r[0].y * v.y + m.r[0].z * v.z, - m.r[1].x * v.x + m.r[1].y * v.y + m.r[1].z * v.z, - m.r[2].x * v.x + m.r[2].y * v.y + m.r[2].z * v.z); -} - -// multiply row vector by a 4x4 matrix -template Vector3 operator*(const Vector3& v, const Matrix4& m) -{ - return Vector3(m.r[0].x * v.x + m.r[1].x * v.y + m.r[2].x * v.z, - m.r[0].y * v.x + m.r[1].y * v.y + m.r[2].y * v.z, - m.r[0].z * v.x + m.r[1].z * v.y + m.r[2].z * v.z); -} - -// multiply column point by a 4x4 matrix; no projection is performed -template Point3 operator*(const Matrix4& m, const Point3& p) -{ - return Point3(m.r[0].x * p.x + m.r[0].y * p.y + m.r[0].z * p.z + m.r[0].w, - m.r[1].x * p.x + m.r[1].y * p.y + m.r[1].z * p.z + m.r[1].w, - m.r[2].x * p.x + m.r[2].y * p.y + m.r[2].z * p.z + m.r[2].w); -} - -// multiply row point by a 4x4 matrix; no projection is performed -template Point3 operator*(const Point3& p, const Matrix4& m) -{ - return Point3(m.r[0].x * p.x + m.r[1].x * p.y + m.r[2].x * p.z + m.r[3].x, - m.r[0].y * p.x + m.r[1].y * p.y + m.r[2].y * p.z + m.r[3].y, - m.r[0].z * p.x + m.r[1].z * p.y + m.r[2].z * p.z + m.r[3].z); -} - -// multiply column vector by a 4x4 matrix -template Vector4 operator*(const Matrix4& m, const Vector4& v) -{ - return Vector4(m.r[0].x * v.x + m.r[0].y * v.y + m.r[0].z * v.z + m.r[0].w * v.w, - m.r[1].x * v.x + m.r[1].y * v.y + m.r[1].z * v.z + m.r[1].w * v.w, - m.r[2].x * v.x + m.r[2].y * v.y + m.r[2].z * v.z + m.r[2].w * v.w, - m.r[3].x * v.x + m.r[3].y * v.y + m.r[3].z * v.z + m.r[3].w * v.w); -} - -// multiply row vector by a 4x4 matrix -template Vector4 operator*(const Vector4& v, const Matrix4& m) -{ - return Vector4(m.r[0].x * v.x + m.r[1].x * v.y + m.r[2].x * v.z + m.r[3].x * v.w, - m.r[0].y * v.x + m.r[1].y * v.y + m.r[2].y * v.z + m.r[3].y * v.w, - m.r[0].z * v.x + m.r[1].z * v.y + m.r[2].z * v.z + m.r[3].z * v.w, - m.r[0].w * v.x + m.r[1].w * v.y + m.r[2].w * v.z + m.r[3].w * v.w); -} - - - -template Matrix4 Matrix4::transpose() const -{ - return Matrix4(Vector4(r[0].x, r[1].x, r[2].x, r[3].x), - Vector4(r[0].y, r[1].y, r[2].y, r[3].y), - Vector4(r[0].z, r[1].z, r[2].z, r[3].z), - Vector4(r[0].w, r[1].w, r[2].w, r[3].w)); -} - - -template Matrix4 operator*(const Matrix4& a, - const Matrix4& b) -{ -#define MATMUL(R, C) (a[R].x * b[0].C + a[R].y * b[1].C + a[R].z * b[2].C + a[R].w * b[3].C) - return Matrix4(Vector4(MATMUL(0, x), MATMUL(0, y), MATMUL(0, z), MATMUL(0, w)), - Vector4(MATMUL(1, x), MATMUL(1, y), MATMUL(1, z), MATMUL(1, w)), - Vector4(MATMUL(2, x), MATMUL(2, y), MATMUL(2, z), MATMUL(2, w)), - Vector4(MATMUL(3, x), MATMUL(3, y), MATMUL(3, z), MATMUL(3, w))); -#undef MATMUL -} - - -template Matrix4 operator+(const Matrix4& a, const Matrix4& b) -{ - return Matrix4(a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]); -} - - -// Compute inverse using Gauss-Jordan elimination; caller is responsible -// for ensuring that the matrix isn't singular. -template Matrix4 Matrix4::inverse() const -{ - Matrix4 a(*this); - Matrix4 b(Matrix4::identity()); - int i, j; - int p; - - for (j = 0; j < 4; j++) - { - p = j; - for (i = j + 1; i < 4; i++) - { - if (fabs(a.r[i][j]) > fabs(a.r[p][j])) - p = i; - } - - // Swap rows p and j - Vector4 t = a.r[p]; - a.r[p] = a.r[j]; - a.r[j] = t; - - t = b.r[p]; - b.r[p] = b.r[j]; - b.r[j] = t; - - T s = a.r[j][j]; // if s == 0, the matrix is singular - a.r[j] *= (1.0f / s); - b.r[j] *= (1.0f / s); - - // Eliminate off-diagonal elements - for (i = 0; i < 4; i++) - { - if (i != j) - { - b.r[i] -= a.r[i][j] * b.r[j]; - a.r[i] -= a.r[i][j] * a.r[j]; - } - } - } - - return b; -} - - -#endif // _VECMATH_H_ diff --git a/src/tools/cmod/cmodsphere/cmodsphere.cpp b/src/tools/cmod/cmodsphere/cmodsphere.cpp index 37f67293d..de22cb173 100644 --- a/src/tools/cmod/cmodsphere/cmodsphere.cpp +++ b/src/tools/cmod/cmodsphere/cmodsphere.cpp @@ -8,7 +8,8 @@ #endif #include -#include +#include + using namespace std; @@ -103,8 +104,8 @@ inline float sampleBilinear(const float samples[], // subdiv is the number of rows in the triangle void triangleSection(unsigned int subdiv, - Vec3f v0, Vec3f v1, Vec3f v2, - Vec2f tex0, Vec2f tex1, Vec2f tex2) + Eigen::Vector3f v0, Eigen::Vector3f v1, Eigen::Vector3f v2, + Eigen::Vector2f tex0, Eigen::Vector2f tex1, Eigen::Vector2f tex2) { float ssamp = (float) (longSamples - 1) + 0.99f; float tsamp = (float) (latSamples - 1) + 0.99f; @@ -116,19 +117,19 @@ void triangleSection(unsigned int subdiv, float u = (i == 0) ? 0.0f : (float) j / (float) i; float v = (float) i / (float) subdiv; - Vec3f w0 = (1.0f - v) * v0 + v * v1; - Vec3f w1 = (1.0f - v) * v0 + v * v2; - Vec3f w = (1.0f - u) * w0 + u * w1; + Eigen::Vector3f w0 = (1.0f - v) * v0 + v * v1; + Eigen::Vector3f w1 = (1.0f - v) * v0 + v * v2; + Eigen::Vector3f w = (1.0f - u) * w0 + u * w1; - Vec2f t((1.0f - u) * tex1.x + u * tex2.x, - (1.0f - v) * tex0.y + v * tex1.y); + Eigen::Vector2f t((1.0f - u) * tex1.x() + u * tex2.x(), + (1.0f - v) * tex0.y() + v * tex1.y()); w.normalize(); if (samples != nullptr) { - float theta = (float) acos(w.y); - float phi = (float) atan2(-w.z, w.x); + float theta = (float) acos(w.y()); + float phi = (float) atan2(-w.z(), w.x()); float s = phi / (2.0f * pi) + 0.5f; float t = theta / pi; @@ -137,8 +138,8 @@ void triangleSection(unsigned int subdiv, w = w * r; } - cout << w.x << " " << w.y << " " << w.z << " " - << t.x << " " << t.y << "\n"; + cout << w.x() << " " << w.y() << " " << w.z() << " " + << t.x() << " " << t.y() << "\n"; } } } @@ -237,62 +238,62 @@ int main(int argc, char* argv[]) cout << "vertices " << vertexCount << "\n"; triangleSection(subdiv, - Vec3f(0.0f, 1.0f, 0.0f), - Vec3f(1.0f, 0.0f, 0.0f), - Vec3f(0.0f, 0.0f, -1.0f), - Vec2f(0.0f, 0.0f), - Vec2f(0.00f, 0.5f), - Vec2f(0.25f, 0.5f)); + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f), + Eigen::Vector2f(0.0f, 0.0f), + Eigen::Vector2f(0.00f, 0.5f), + Eigen::Vector2f(0.25f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, 1.0f, 0.0f), - Vec3f(0.0f, 0.0f, 1.0f), - Vec3f(1.0f, 0.0f, 0.0f), - Vec2f(0.0f, 0.0f), - Vec2f(0.75f, 0.5f), - Vec2f(1.00f, 0.5f));; + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector2f(0.0f, 0.0f), + Eigen::Vector2f(0.75f, 0.5f), + Eigen::Vector2f(1.00f, 0.5f));; triangleSection(subdiv, - Vec3f(0.0f, 1.0f, 0.0f), - Vec3f(-1.0f, 0.0f, 0.0f), - Vec3f(0.0f, 0.0f, 1.0f), - Vec2f(0.0f, 0.0f), - Vec2f(0.50f, 0.5f), - Vec2f(0.75f, 0.5f)); + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector2f(0.0f, 0.0f), + Eigen::Vector2f(0.50f, 0.5f), + Eigen::Vector2f(0.75f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, 1.0f, 0.0f), - Vec3f(0.0f, 0.0f, -1.0f), - Vec3f(-1.0f, 0.0f, 0.0f), - Vec2f(0.0f, 0.0f), - Vec2f(0.25f, 0.5f), - Vec2f(0.50f, 0.5f)); + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector2f(0.0f, 0.0f), + Eigen::Vector2f(0.25f, 0.5f), + Eigen::Vector2f(0.50f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, -1.0f, 0.0f), - Vec3f(0.0f, 0.0f, -1.0f), - Vec3f(1.0f, 0.0f, 0.0f), - Vec2f(0.0f, 1.0f), - Vec2f(0.25f, 0.5f), - Vec2f(0.00f, 0.5f)); + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f), + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector2f(0.0f, 1.0f), + Eigen::Vector2f(0.25f, 0.5f), + Eigen::Vector2f(0.00f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, -1.0f, 0.0f), - Vec3f(1.0f, 0.0f, 0.0f), - Vec3f(0.0f, 0.0f, 1.0f), - Vec2f(0.0f, 1.0f), - Vec2f(1.00f, 0.5f), - Vec2f(0.75f, 0.5f)); + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector2f(0.0f, 1.0f), + Eigen::Vector2f(1.00f, 0.5f), + Eigen::Vector2f(0.75f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, -1.0f, 0.0f), - Vec3f(0.0f, 0.0f, 1.0f), - Vec3f(-1.0f, 0.0f, 0.0f), - Vec2f(0.0f, 1.0f), - Vec2f(0.75f, 0.5f), - Vec2f(0.50f, 0.5f)); + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector2f(0.0f, 1.0f), + Eigen::Vector2f(0.75f, 0.5f), + Eigen::Vector2f(0.50f, 0.5f)); triangleSection(subdiv, - Vec3f(0.0f, -1.0f, 0.0f), - Vec3f(-1.0f, 0.0f, 0.0f), - Vec3f(0.0f, 0.0f, -1.0f), - Vec2f(0.0f, 1.0f), - Vec2f(0.50f, 0.5f), - Vec2f(0.25f, 0.5f)); + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f), + Eigen::Vector2f(0.0f, 1.0f), + Eigen::Vector2f(0.50f, 0.5f), + Eigen::Vector2f(0.25f, 0.5f)); cout << "trilist 0 " << triangleCount * 3 << "\n";