📄 matrix.h
字号:
// For orthogonal matrices, I belive this also gives you the inverse.
void Transpose()
{
float f;
SWAP(f12, f21, f);
SWAP(f13, f31, f);
SWAP(f14, f41, f);
SWAP(f23, f32, f);
SWAP(f24, f42, f);
SWAP(f34, f43, f);
}
};
/*******************************************************************************
* Class: CLine
********************************************************************************
* This class implements a 3-dimensional line. It is initialized with two points
* on the line. It can be used to find the shortest distance between two lines or
* between a point and a line, which can be useful for things like collision
* detection.
*******************************************************************************/
class CLine
{
public:
CVector m_vStart;
CVector m_vDir;
CLine() {}
CLine(CVector &p1, CVector &p2) { Init(p1, p2); }
void Init(CVector &p1, CVector &p2)
{
m_vStart = p1;
m_vDir = p2 - p1;
m_vDir.Normalize();
}
CVector ClosestPoint(CVector &p) { return m_vStart + (m_vDir * ((m_vDir | p) - (m_vDir | m_vStart))); }
float Distance(CVector &p) { return p.Distance(ClosestPoint(p)); }
float Distance(CLine &l)
{
CVector v = m_vDir ^ l.m_vDir;
v.Normalize();
return Abs((m_vStart - l.m_vStart) | v);
}
/*
// Calculate the shortest distance between a point to a line
(known) point p = p, or (px, py, pz)
(known) point on line = a, or (ax, ay, az)
(known) direction vector of line = n, or (nx, ny, nz)
(unknown) point on line closest to p = P, or (Px, Py, Pz)
// The vector from p to P must be perpendicular to n, which gives us:
(px-Px, py-Py, pz-Pz) . (nx, ny, nz) = 0
// To make sure P is actually on the line, it must also fit this line equation:
(ax, ay, az) + k * (nx, ny, nz) = (Px, Py, Pz)
// Expanding those equations gives us this:
nx * (px - Px) + ny * (py - Py) + nz * (pz - Pz) = 0
// and these:
Px = ax + k * nx
Py = ay + k * ny
Pz = az + k * nz
// Substituting and solving for k gives us:
nx * (px - (ax + k * nx)) + ny * (py - (ay + k * ny)) + nz * (pz - (az + k * nz)) = 0
nx*px - nx*ax - k*nx2 + ny*py - ny*ay - k*ny2 + nz*pz - nz*az - k*nz2 = 0
k * (nx2 + ny2 + nz2) = nx*px - nx*ax + ny*py - ny*ay + nz*pz - nz*az
k = (nx*px + ny*py + nz*pz - nx*ax - ny*ay - nz*az) / (nx2 + ny2 + nz2)
k = (n.p - n.a) / (n.n)
// If n is normalized then n.n will always equal 1, so...
P = a + (n.p - n.a) * n
distance = | P - p |
*/
/*
// Calculate the shortest distance between two lines (call them line A and line B)
(known) point on A = AP, or (APx, APy, APz)
(known) direction vector of A = AN, or (ANx, ANy, ANz)
(known) point on B = BP, or (BPx, BPy, BPz)
(known) direction vector of B = BN, or (BNx, BNy, BNz)
(unknown) point on A closest to B = P1, or (P1x, P1y, P1z)
(unknown) point on B closest to A = P2, or (P2x, P2y, P2z)
// The vector from P1 to P2 is perpendicular to both lines:
(P2x - P1x, P2y - P1y, P2z - P1z) . (ANx, ANy, ANz) = 0
(P2x - P1x, P2y - P1y, P2z - P1z) . (BNx, BNy, BNz) = 0
// To find P1 and P2, we need to find values k and l for the line equations:
(APx, APy, APz) + k * (ANx, ANy, ANz) = (P1x, P1y, P1z)
(BPx, BPy, BPz) + l * (BNx, BNy, BNz) = (P2x, P2y, P2z)
// Expanding these equations gives us these:
ANx * (P2x - P1x) + ANy * (P2y - P1y) + ANz * (P2z - P1z) = 0
BNx * (P2x - P1x) + BNy * (P2y - P1y) + BNz * (P2z - P1z) = 0
// And these:
P1x = APx + k * ANx
P1y = APy + k * ANy
P1z = APz + k * ANz
P2x = BPx + l * BNx
P2y = BPy + l * BNy
P2z = BPz + l * BNz
// Substituting and solving for k and l gives us:
ANx * ((BPx + l * BNx) - (APx + k * ANx)) + ANy * ((BPy + l * BNy) - (APy + k * ANy)) + ANz * ((BPz + l * BNz) - (APz + k * ANz)) = 0
BNx * ((BPx + l * BNx) - (APx + k * ANx)) + BNy * ((BPy + l * BNy) - (APy + k * ANy)) + BNz * ((BPz + l * BNz) - (APz + k * ANz)) = 0
(ANx*BPx + l*ANx*BNx - ANx*APx - k*ANx2) + (ANy*BPy + l*ANy*BNy - ANy*APy - k*ANy2) + (ANz*BPz + l*ANz*BNz - ANz*APz - k*ANz2) = 0
(BNx*BPx + l*BNx2 - BNx*APx - k*ANx*BNx) + (BNy*BPy + l*BNy2 - BNy*APy - k*ANy*BNy) + (BNz*BPz + l*BNz2 - BNz*APz - k*ANz*BNz) = 0
k = (ANx*BPx + l*ANx*BNx - ANx*APx + ANy*BPy + l*ANy*BNy - ANy*APy + ANz*BPz + l*ANz*BNz - ANz*APz) / (ANx2 + ANy2 + ANz2)
k = (BNx*BPx + l*BNx2 - BNx*APx + BNy*BPy + l*BNy2 - BNy*APy + BNz*BPz + l*BNz2 - BNz*APz) / (ANx*BNx + ANy*BNy + ANz*BNz)
(ANx*BNx + ANy*BNy + ANz*BNz) * (ANx*BPx + l*ANx*BNx - ANx*APx + ANy*BPy + l*ANy*BNy - ANy*APy + ANz*BPz + l*ANz*BNz - ANz*APz) = (ANx2 + ANy2 + ANz2) * (BNx*BPx + l*BNx2 - BNx*APx + BNy*BPy + l*BNy2 - BNy*APy + BNz*BPz + l*BNz2 - BNz*APz)
(ANx*BNx*ANy*BPy + ANx*BNx*l*ANy*BNy - ANx*BNx*ANy*APy + ANx*BNx*ANz*BPz + ANx*BNx*l*ANz*BNz - ANx*BNx*ANz*APz) +
(ANy*BNy*ANx*BPx + ANy*BNy*l*ANx*BNx - ANy*BNy*ANx*APx + ANy*BNy*ANz*BPz + ANy*BNy*l*ANz*BNz - ANy*BNy*ANz*APz) +
(ANz*BNz*ANx*BPx + ANz*BNz*l*ANx*BNx - ANz*BNz*ANx*APx + ANz*BNz*ANy*BPy + ANz*BNz*l*ANy*BNy - ANz*BNz*ANy*APy) =
(ANx2*BNy*BPy + ANx2*l*BNy2 - ANx2*BNy*APy + ANx2*BNz*BPz + ANx2*l*BNz2 - ANx2*BNz*APz) +
(ANy2*BNx*BPx + ANy2*l*BNx2 - ANy2*BNx*APx + ANy2*BNz*BPz + ANy2*l*BNz2 - ANy2*BNz*APz) +
(ANz2*BNx*BPx + ANz2*l*BNx2 - ANz2*BNx*APx + ANz2*BNy*BPy + ANz2*l*BNy2 - ANz2*BNy*APy)
l * (2*ANx*BNx*ANy*BNy + 2*ANz*BNz*ANx*BNx + 2*ANy*BNy*ANz*BNz - ANx2*BNy2 - ANy2*BNx2 - ANz2*BNx2 - ANx2*BNz2 - ANy2*BNz2 - ANz2*BNy2) =
(ANx2*BNy*BPy - ANx2*BNy*APy + ANx2*BNz*BPz - ANx2*BNz*APz) +
(ANy2*BNx*BPx - ANy2*BNx*APx + ANy2*BNz*BPz - ANy2*BNz*APz) +
(ANz2*BNx*BPx - ANz2*BNx*APx + ANz2*BNy*BPy - ANz2*BNy*APy) +
(-ANx*BNx*ANy*BPy + ANx*BNx*ANy*APy - ANx*BNx*ANz*BPz + ANx*BNx*ANz*APz) +
(-ANy*BNy*ANx*BPx + ANy*BNy*ANx*APx - ANy*BNy*ANz*BPz + ANy*BNy*ANz*APz) +
(-ANz*BNz*ANx*BPx + ANz*BNz*ANx*APx - ANz*BNz*ANy*BPy + ANz*BNz*ANy*APy)
// Next think about using AN ^ BN to get the direction vector between P1 and P2
// Normalize that vector and call it n
// The distance between the two lines might be the absolute value of (P2-P1).n
// That still doesn't tell us where the closest points on the line are, but it might help find them
*/
};
/*******************************************************************************
* Class: CPlane
********************************************************************************
* This class implements a 3-dimensional plane equation. It can be initialized
* with 3 points on the plane or a normal and a point on the plane, and it can
* be used to find the distance from a point to a plane. I plan to add more to
* this class later to make it more useful.
*******************************************************************************/
class CPlane
{
public:
CVector m_vNormal;
float D;
CPlane() {}
void Init(CVector &p1, CVector &p2, CVector &p3)
{ // Initializes the plane based on three points in the plane
Init(NormalVector(p1, p2, p3), p1);
}
void Init(CVector &vNormal, CVector &p)
{ // Initializes the plane based on a normal and a point in the plane
m_vNormal = vNormal;
D = -(p | m_vNormal);
}
bool operator==(CPlane &p) { return (Abs(D - p.D) <= DELTA && m_vNormal == p.m_vNormal); }
bool operator!=(CPlane &p) { return !(*this == p); }
float Distance(CVector &p)
{ // Returns the distance between the plane and point p
return (m_vNormal | p) + D; // A positive, 0, or negative result indicates the point is in front of, on, or behind the plane
}
bool Intersection(CVector &vPos, CVector &vDir)
{ // Returns true if the line intersects the plane and changes vPos to the location of the intersection
float f = m_vNormal | vDir;
if(ABS(f) < DELTA)
return false;
vPos -= vDir * (Distance(vPos) / f);
return true;
}
};
/****************************************************************************
* Class: C3DObject
*****************************************************************************
* This class represents a basic 3D object in the scene. It has a 3D position,
* orientation, velocity, and a parent which provides its frame of reference
* in the scene.
* Note: This class is derived from CQuaternion so it will inherit useful
* functions like Rotate(), GetViewAxis(), GetUpAxis(), and GetRightAxis().
****************************************************************************/
#define MIN_DISTANCE 0.01
#define MAX_DISTANCE 1000.0 // Distance to desired far clipping plane
#define MAX_DISCERNABLE 1000000.0 // Beyond this distance, everything is rendered at MAX_DISTANCE
#define HALF_MAX (MAX_DISTANCE*0.5) // Everything between HALF_MAX and MAX_DISCERNABLE is scaled exponentially between HALF_MAX and MAX_DISTANCE
class C3DObject : public CQuaternion
{
protected:
C3DObject *m_pParent; // The object's parent
float m_fMass; // The object's mass (kg)
CDoubleVector m_vPosition; // The object's position (km)
CVector m_vVelocity; // The object's velocity (km/s)
public:
C3DObject() : CQuaternion(0.0f, 0.0f, 0.0f, 1.0f), m_vPosition(0.0f), m_vVelocity(0.0f)
{
m_fMass = 0.0f;
}
void SetParent(C3DObject *p) { m_pParent = p; }
C3DObject *GetParent() { return m_pParent; }
void SetMass(float f) { m_fMass = f; }
float GetMass() { return m_fMass; }
void SetPosition(CDoubleVector &v) { m_vPosition = v; }
CDoubleVector GetPosition() { return m_vPosition; }
void SetVelocity(CVector &v) { m_vVelocity = v; }
CVector GetVelocity() { return m_vVelocity; }
CMatrix GetViewMatrix()
{
// Don't use the normal view matrix because it causes precision problems if the camera is too far away from the origin.
// Instead, pretend the camera is at the origin and offset all model matrices by subtracting the camera's position.
CMatrix m = *this;
m.Transpose();
return m;
}
CMatrix GetModelMatrix(C3DObject *pCamera)
{
// Don't use the normal model matrix because it causes precision problems if the camera and model are too far away from the origin.
// Instead, pretend the camera is at the origin and offset all model matrices by subtracting the camera's position.
CMatrix m;
m.ModelMatrix(*this, m_vPosition-pCamera->m_vPosition);
return m;
}
float ScaleModelMatrix(CMatrix &m, double dDistance, double dFactor)
{
// This code scales the object's size and distance to the camera down when it's too far away.
// This solves a problem with many video card drivers where objects too far away aren't rendering properly.
// It also alleviates the Z-buffer precision problem caused by having your near and far clipping planes too far apart.
if(dDistance > HALF_MAX)
{
dFactor *= (dDistance >= MAX_DISCERNABLE) ? MAX_DISTANCE : HALF_MAX + HALF_MAX * (1.0 - exp(-2.5 * dDistance / MAX_DISCERNABLE));
dFactor /= dDistance;
}
float fFactor = (float)dFactor;
if(fFactor <= 1.0f)
{
m.f41 *= fFactor;
m.f42 *= fFactor;
m.f43 *= fFactor;
m.Scale(fFactor, fFactor, fFactor);
}
return fFactor;
}
float GetScaledModelMatrix(CMatrix &m, C3DObject *pCamera, double dFactor=1.0)
{
CDoubleVector vPos = m_vPosition - pCamera->m_vPosition;
double dDistance = vPos.Magnitude();
m.ModelMatrix(*this, vPos);
return ScaleModelMatrix(m, dDistance, dFactor);
}
float GetScaledBillboardMatrix(CMatrix &m, C3DObject *pCamera, float fSize, double dFactor=1.0)
{
CDoubleVector vPos = m_vPosition - pCamera->m_vPosition;
double dDistance = vPos.Magnitude();
CVector vView = vPos / dDistance;
CVector vRight = vView ^ CVector(0, 1, 0);
vRight.Normalize();
CVector vUp = vRight ^ vView;
vUp.Normalize();
m.ModelMatrix(vPos, vView, vUp, vRight);
float fFactor = ScaleModelMatrix(m, dDistance, dFactor);
m.Scale(fSize, fSize, fSize);
return fFactor;
}
float DistanceTo(C3DObject &obj) { return (float)m_vPosition.Distance(obj.m_vPosition); }
CVector VectorTo(C3DObject &obj) { return obj.m_vPosition - m_vPosition; }
float GravityPull(float fDistSquared) { return (GRAVCONST * m_fMass) / fDistSquared; }
CVector GravityVector(C3DObject &obj)
{
CVector vGravity = obj.VectorTo(*this);
float fDistSquared = vGravity.MagnitudeSquared();
vGravity *= GravityPull(fDistSquared) / sqrtf(fDistSquared);
return vGravity;
}
void Accelerate(CVector &vAccel, float fSeconds, float fResistance=0)
{
m_vVelocity += vAccel * fSeconds;
if(fResistance > DELTA)
m_vVelocity *= 1.0f - fResistance * fSeconds;
m_vPosition += m_vVelocity * fSeconds;
}
// Kinetic energy (Joules) = 1/2 * mass * velocity^2 (mass in kg, velocity in m/s)
// Kinetic energy (kilotons) = (KE in Joules) / (4.185e12 joules/KT)
float GetKEJoules(C3DObject &obj) { return 500000.0f * m_fMass * (m_vVelocity-obj.m_vVelocity).MagnitudeSquared(); }
float GetKEKilotons(C3DObject &obj) { return 1.19474e-7f * m_fMass * (m_vVelocity-obj.m_vVelocity).MagnitudeSquared(); }
};
/*
class CPlanetoid : public C3DObject
{
protected:
float m_fRadius;
float m_fDensity;
float m_fImpactFactor;
public:
CPlanetoid() : C3DObject() {}
void Init(float fMass, float fRadius)
{
m_fRadius = fRadius;
m_fMass = fMass;
float f = 1000*m_fRadius; // Get radius in meters and its square
float fDensity = (0.75f*INV_PI)*m_fMass/(f*f*f);// Get density (kg/m^3)
float fGravity = 1000*GravityPull(m_fRadius*m_fRadius); // Get pull of gravity at surface (m/s^2)
m_fImpactFactor = 0.037f * powf(9.8f/fGravity, 0.166667f) * powf(1.8e3f/fDensity, 0.294118f);
}
// Impact crater radius (km) = m_fImpactFactor * (KE in kilotons)^0.294118
float GetCraterRadius(C3DObject &obj)
{
return m_fImpactFactor * powf(GetKEKilotons(obj), 0.294118f);
}
};
*/
#endif //__Matrix_h__
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -