OU/CWE-ou-minkata
1
0
Fork 0
mirror of https://foundry.openuru.org/gitblit/r/CWE-ou-minkata.git synced 2025-07-17 10:52:46 +00:00
Code Releases Activity

Initial Commit of CyanWorlds.com Engine Open Source Client/Plugin

Browse Source
This commit is contained in:
JWPlatt
2011-03-12 12:34:52 -05:00
commit a20a222fc2
3976 changed files with 1301355 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

426
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/hsAffineParts.cpp Normal file
View File

@ -0,0 +1,426 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
#include "HeadSpin.h"
#include "hsAffineParts.h"
#include "../plInterp/hsInterp.h"
#include "hsStream.h"
#include "plProfile.h"
#define PL_OPTIMIZE_COMPOSE
inline void QuatTo3Vectors(const hsQuat& q, hsVector3* const v)
{
v[0][0] = 1.0f - 2.0f*q.fY*q.fY - 2.0f*q.fZ*q.fZ;
v[0][1] = 2.0f*q.fX*q.fY - 2.0f*q.fW*q.fZ;
v[0][2] = 2.0f*q.fX*q.fZ + 2.0f*q.fW*q.fY;
v[1][0] = 2.0f*q.fX*q.fY + 2.0f*q.fW*q.fZ;
v[1][1] = 1.0f - 2.0f*q.fX*q.fX - 2.0f*q.fZ*q.fZ;
v[1][2] = 2.0f*q.fY*q.fZ - 2.0f*q.fW*q.fX;
v[2][0] = 2.0f*q.fX*q.fZ - 2.0f*q.fW*q.fY;
v[2][1] = 2.0f*q.fY*q.fZ + 2.0f*q.fW*q.fX;
v[2][2] = 1.0f - 2.0f*q.fX*q.fX - 2.0f*q.fY*q.fY;
}
inline void QuatTo3VectorsTranspose(const hsQuat& q, hsVector3* const v)
{
v[0][0] = 1.0f - 2.0f*q.fY*q.fY - 2.0f*q.fZ*q.fZ;
v[1][0] = 2.0f*q.fX*q.fY - 2.0f*q.fW*q.fZ;
v[2][0] = 2.0f*q.fX*q.fZ + 2.0f*q.fW*q.fY;
v[0][1] = 2.0f*q.fX*q.fY + 2.0f*q.fW*q.fZ;
v[1][1] = 1.0f - 2.0f*q.fX*q.fX - 2.0f*q.fZ*q.fZ;
v[2][1] = 2.0f*q.fY*q.fZ - 2.0f*q.fW*q.fX;
v[0][2] = 2.0f*q.fX*q.fZ - 2.0f*q.fW*q.fY;
v[1][2] = 2.0f*q.fY*q.fZ + 2.0f*q.fW*q.fX;
v[2][2] = 1.0f - 2.0f*q.fX*q.fX - 2.0f*q.fY*q.fY;
}
//
// Constructors
// Convert from Gems struct for now
//
hsAffineParts::hsAffineParts(gemAffineParts *ap)
{
AP_SET((*this), (*ap));
}
//
//
//
hsAffineParts::hsAffineParts()
{
}
//
//
//
void hsAffineParts::Reset()
{
fT.Set(0,0,0);
fQ.Identity();
fU.Identity();
fK.Set(1,1,1);
fF = 1.0;
}
plProfile_CreateTimer("Compose", "Affine", Compose);
plProfile_CreateTimer("ComposeInv", "Affine", ComposeInv);
//
// Create an affine matrix from the various parts
//
// AffineParts:
// Vector t; /* Translation components */
// Quat q; /* Essential rotation */
// Quat u; /* Stretch rotation */
// Vector k; /* Stretch factors */
// float f; /* Sign of determinant */
//
// A matrix M is decomposed by : M = T F R U K Utranspose.
// T is the translate mat.
// F is +-Identity (to flip the rotation or not).
// R is the rot matrix.
// U is the stretch matrix.
// K is the scale factor matrix.
//
void hsAffineParts::ComposeMatrix(hsMatrix44 *out) const
{
plProfile_BeginTiming(Compose);
#ifndef PL_OPTIMIZE_COMPOSE
// Built U matrix
hsMatrix44 U;
fU.MakeMatrix(&U);
// Build scale factor matrix
hsMatrix44 K;
K.MakeScaleMat(&fK);
// Build Utranspose matrix
hsMatrix44 Utp;
U.GetTranspose(&Utp);
// Build R matrix
hsMatrix44 R;
fQ.MakeMatrix(&R);
// Build flip matrix
// hsAssert(fF == 1.0 || fF == -1.0, "Invalid flip portion of affine parts");
hsMatrix44 F;
if (fF==-1.0)
{
hsVector3 s;
s.Set(-1,-1,-1);
F.MakeScaleMat(&s);
}
else
F.Reset();
// Build translate matrix
hsMatrix44 T;
T.MakeTranslateMat(&fT);
//
// Concat mats
//
*out = K * Utp;
*out = U * (*out);
*out = R * (*out); // Q
*out = F * (*out);
*out = T * (*out); // Translate happens last
#else // PL_OPTIMIZE_COMPOSE
// M = T F R U K Ut,
// but these are mostly very sparse matrices. So rather
// than construct the full 6 matrices and concatenate them,
// we'll work out by hand what the non-zero results will be.
// T = |1 0 0 Tx|
// |0 1 0 Ty|
// |0 0 1 Tz|
// F = |f 0 0 0|
// |0 f 0 0|
// |0 0 f 0|, where f is either 1 or -1
// R = |R00 R01 R02 0|
// |R10 R11 R12 0|
// |R20 R21 R22 0|
// U = |U00 U01 U02 0|
// |U10 U11 U12 0|
// |U20 U21 U22 0|
// K = |Sx 0 0 0|
// |0 Sy 0 0|
// |0 0 Sz 0|
// Ut = |U00 U10 U20 0|
// |U01 U11 U21 0|
// |U02 U12 U22 0|, where Uij is from matrix U
//
// So, K * Ut =
// |Sx*U00 Sx*U10 Sx*U20 0|
// |Sy*U01 Sy*U11 Sy*U21 0|
// |Sz*U02 Sz*U12 Sz*U22 0|
//
// U * (K * Ut) =
// | U0 dot S*U0 U0 dot S*U1 U0 dot S*U2 0|
// | U1 dot S*U0 U1 dot S*U1 U1 dot S*U2 0|
// | U2 dot S*U0 U2 dot S*U1 U2 dot S*U2 0|
//
// Let's call that matrix UK
//
// Now R * U * K * Ut = R * UK =
// | R0 dot UKc0 R0 dot UKc1 R0 dot UKc2 0|
// | R1 dot UKc0 R1 dot UKc1 R1 dot UKc2 0|
// | R2 dot UKc0 R2 dot UKc1 R2 dot UKc2 0|, where UKci is column i from UK
//
// if f is -1, we negate the matrix we have so far, else we don't. We can
// accomplish this cleanly by just negating the scale vector S if f == -1.
//
// Since the translate is last, we can just stuff it into the 4th column.
//
// Since we only ever use UK as column vectors, we'll just construct it
// into 3 vectors representing the columns.
//
// The quat MakeMatrix function is pretty efficient, but it does a little more work
// than it has to filling out the whole matrix when we only need the 3x3 rotation,
// and we'd rather have it in the form of vectors anyway, so we'll use our own
// quat to 3 vectors function here.
hsVector3 U[3];
QuatTo3Vectors(fU, U);
int i, j;
hsVector3 UKt[3];
for( i = 0; i < 3; i++ )
{
for( j = 0; j < 3; j++ )
{
// SU[j] = (fK.fX * U[j].fX, fK.fY * U[j].fY, fK.fZ * U[j].fZ)
UKt[j][i] = U[i].fX * fK.fX * U[j].fX
+ U[i].fY * fK.fY * U[j].fY
+ U[i].fZ * fK.fZ * U[j].fZ;
}
}
hsVector3 R[3];
QuatTo3Vectors(fQ, R);
hsScalar f = fF < 0 ? -1.f : 1.f;
for( i = 0; i < 3; i++ )
{
for( j = 0; j < 3; j++ )
{
out->fMap[i][j] = R[i].InnerProduct(UKt[j]) * f;
}
out->fMap[i][3] = fT[i];
}
out->fMap[3][0] = out->fMap[3][1] = out->fMap[3][2] = 0.f;
out->fMap[3][3] = 1.f;
out->NotIdentity();
#endif // PL_OPTIMIZE_COMPOSE
plProfile_EndTiming(Compose);
}
void hsAffineParts::ComposeInverseMatrix(hsMatrix44 *out) const
{
plProfile_BeginTiming(Compose);
#ifndef PL_OPTIMIZE_COMPOSE
// Built U matrix
hsMatrix44 U;
fU.Conjugate().MakeMatrix(&U);
// Build scale factor matrix
hsMatrix44 K;
hsVector3 invK;
invK.Set(hsScalarInvert(fK.fX),hsScalarInvert(fK.fY),hsScalarInvert(fK.fZ));
K.MakeScaleMat(&invK);
// Build Utranspose matrix
hsMatrix44 Utp;
U.GetTranspose(&Utp);
// Build R matrix
hsMatrix44 R;
fQ.Conjugate().MakeMatrix(&R);
// Build flip matrix
// hsAssert(fF == 1.0 || fF == -1.0, "Invalid flip portion of affine parts");
hsMatrix44 F;
if (fF==-1.0)
{
hsVector3 s;
s.Set(-1,-1,-1);
F.MakeScaleMat(&s);
}
else
F.Reset();
// Build translate matrix
hsMatrix44 T;
T.MakeTranslateMat(&-fT);
//
// Concat mats
//
*out = Utp * K;
*out = (*out) * U;
*out = (*out) * R;
*out = (*out) * F;
*out = (*out) * T;
#else // PL_OPTIMIZE_COMPOSE
// Same kind of thing here, except now M = Ut K U R F T
// and again
// T = |1 0 0 Tx|
// |0 1 0 Ty|
// |0 0 1 Tz|
// F = |f 0 0 0|
// |0 f 0 0|
// |0 0 f 0|, where f is either 1 or -1
// R = |R00 R01 R02 0|
// |R10 R11 R12 0|
// |R20 R21 R22 0|
// U = |U00 U01 U02 0|
// |U10 U11 U12 0|
// |U20 U21 U22 0|
// K = |Sx 0 0 0|
// |0 Sy 0 0|
// |0 0 Sz 0|
// Ut = |U00 U10 U20 0|
// |U01 U11 U21 0|
// |U02 U12 U22 0|, where Uij is from matrix U
//
// So, Ut * K =
// |U00*Sx U10*Sy U20*Sz 0|
// |U01*Sx U11*Sy U21*Sz 0|
// |U02*Sx U12*Sy U22*Sz 0|
//
// (Ut * K) * U = UK =
// |Ut0*S dot Ut0 Ut0*S dot Ut1 Ut0*S dot Ut2 0|
// |Ut1*S dot Ut0 Ut1*S dot Ut1 Ut1*S dot Ut2 0|
// |Ut2*S dot Ut0 Ut2*S dot Ut1 Ut2*S dot Ut2 0|
//
// (((Ut * K) * U) * R)[i][j] = UK[i] dot Rc[j]
//
// Again we'll stuff the flip into the scale.
//
// Now, because the T is on the other end of the concat (closest
// to the vertex), we can't just stuff it in. If Mr is the
// rotation part of the final matrix (Ut * K * U * R * F), then
// the translation components M[i][3] = Mr[i] dot T.
//
//
hsVector3 Ut[3];
QuatTo3VectorsTranspose(fU.Conjugate(), Ut);
int i, j;
hsVector3 invK;
invK.Set(hsScalarInvert(fK.fX),hsScalarInvert(fK.fY),hsScalarInvert(fK.fZ));
hsVector3 UK[3];
for( i = 0; i < 3; i++ )
{
for( j = 0; j < 3; j++ )
{
// SUt[i] = (Ut[i].fX * invK.fX, Ut[i].fY * invK.fY, Ut[i].fZ * invK.fZ)
// So SUt[i].InnerProduct(Ut[j]) ==
// Ut[i].fX * invK.fX * Ut[j].fX
// + Ut[i].fY * invK.fY * Ut[j].fY
// + Ut[i].fZ * invK.fZ * Ut[j].fZ
UK[i][j] = Ut[i].fX * invK.fX * Ut[j].fX
+ Ut[i].fY * invK.fY * Ut[j].fY
+ Ut[i].fZ * invK.fZ * Ut[j].fZ;
}
}
hsVector3 Rt[3];
QuatTo3VectorsTranspose(fQ.Conjugate(), Rt);
hsScalar f = fF < 0 ? -1.f : 1.f;
for( i = 0; i < 3; i++ )
{
for( j = 0; j < 3; j++ )
{
out->fMap[i][j] = UK[i].InnerProduct(Rt[j]) * f;
}
out->fMap[i][3] = -(fT.InnerProduct((hsPoint3*)(&out->fMap[i])));
}
out->fMap[3][0] = out->fMap[3][1] = out->fMap[3][2] = 0.f;
out->fMap[3][3] = 1.f;
out->NotIdentity();
#endif // PL_OPTIMIZE_COMPOSE
plProfile_EndTiming(Compose);
}
//
// Given 2 affineparts structs and a p value (between 0-1),
// compute a new affine parts.
//
void hsAffineParts::SetFromInterp(const hsAffineParts &ap1, const hsAffineParts &ap2, float p)
{
hsAssert(p>=0.0 && p<=1.0, "Interpolate param must be 0-1");
#if 0
// Debug
float rad1,rad2, rad3;
hsVector3 axis1, axis2, axis3;
k1->fQ.GetAngleAxis(&rad1, &axis1);
k2->fQ.GetAngleAxis(&rad2, &axis2);
fQ.GetAngleAxis(&rad3, &axis3);
#endif
hsInterp::LinInterp(&ap1, &ap2, p, this);
}
//
// Read
//
void hsAffineParts::Read(hsStream *stream)
{
fT.Read(stream);
fQ.Read(stream);
fU.Read(stream);
fK.Read(stream);
fF = stream->ReadSwapFloat();
}
//
// Write
//
void hsAffineParts::Write(hsStream *stream)
{
fT.Write(stream);
fQ.Write(stream);
fU.Write(stream);
fK.Write(stream);
stream->WriteSwapFloat(fF);
}

73
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/hsAffineParts.h Normal file
View File

@ -0,0 +1,73 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
#ifndef HSAFFINEPARTS_inc
#define HSAFFINEPARTS_inc
#include "hsGeometry3.h"
#include "hsQuat.h"
#include "mat_decomp.h"
class hsAffineParts
{
public:
// Constructors
hsAffineParts(gemAffineParts *); // Convert from Gems struct for now
hsAffineParts();
void Reset();
hsVector3 fT; /* Translation components */
hsQuat fQ; /* Essential rotation */
hsQuat fU; /* Stretch rotation */
hsVector3 fK; /* Stretch factors */
float fF; /* Sign of determinant */
void ComposeMatrix(hsMatrix44 *out) const;
void ComposeInverseMatrix(hsMatrix44 *out) const;
void SetFromInterp(const hsAffineParts &ap1, const hsAffineParts &ap2, float t);
void Read(hsStream *);
void Write(hsStream *);
int operator==(const hsAffineParts& a) const
{ return (fT == a.fT && fQ == a.fQ && fU == a.fU && fK == a.fK && fF == a.fF); }
};
//
// General set macro can also be used for 3DSMax struct
//
#define AP_SET(dst, src) \
{ \
dst.fT.Set(src.t.x, src.t.y, src.t.z); \
dst.fQ.Set(src.q.x, src.q.y, src.q.z, src.q.w); \
dst.fU.Set(src.u.x, src.u.y, src.u.z, src.u.w); \
dst.fK.Set(src.k.x, src.k.y, src.k.z); \
dst.fF = src.f; \
}
#endif

212
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/hsEuler.cpp Normal file
View File

@ -0,0 +1,212 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
//
//////////////////////////////////////////////////////////////////////////
// EULER STUFF
// See Gems IV, Ken Shoemake
//////////////////////////////////////////////////////////////////////////
//
#include <float.h> // for FLT_EPSILON
#include "hsEuler.h"
#include "hsQuat.h"
#include "hsMatrix44.h"
enum QuatPart
{
X, Y, Z, W
};
//
// Construct quaternion from Euler angles (in radians).
//
void hsEuler::GetQuat(hsQuat* qu)
{
double a[3], ti, tj, th, ci, cj, ch, si, sj, sh, cc, cs, sc, ss;
int i,j,k,h,n,s,f;
hsEuler ea=*this; // copy
EulGetOrd(ea.fOrder,i,j,k,h,n,s,f);
if (f==EulFrmR)
{
hsScalar t = ea.fX; ea.fX = ea.fZ; ea.fZ = t;
}
if (n==EulParOdd)
ea.fY = -ea.fY;
ti = ea.fX*0.5; tj = ea.fY*0.5; th = ea.fZ*0.5;
ci = cos(ti); cj = cos(tj); ch = cos(th);
si = sin(ti); sj = sin(tj); sh = sin(th);
cc = ci*ch; cs = ci*sh; sc = si*ch; ss = si*sh;
if (s==EulRepYes)
{
a[i] = cj*(cs + sc); /* Could speed up with */
a[j] = sj*(cc + ss); /* trig identities. */
a[k] = sj*(cs - sc);
qu->fW = static_cast<float>(cj*(cc - ss));
}
else
{
a[i] = cj*sc - sj*cs;
a[j] = cj*ss + sj*cc;
a[k] = cj*cs - sj*sc;
qu->fW = static_cast<float>(cj*cc + sj*ss);
}
if (n==EulParOdd)
a[j] = -a[j];
qu->fX = static_cast<float>(a[X]);
qu->fY = static_cast<float>(a[Y]);
qu->fZ = static_cast<float>(a[Z]);
}
//
// Construct matrix from Euler angles (in radians).
//
void hsEuler::GetMatrix44(hsMatrix44* mat)
{
double ti, tj, th, ci, cj, ch, si, sj, sh, cc, cs, sc, ss;
int i,j,k,h,n,s,f;
hsEuler ea=*this; // copy
EulGetOrd(ea.fOrder,i,j,k,h,n,s,f);
if (f==EulFrmR)
{
hsScalar t = ea.fX; ea.fX = ea.fZ; ea.fZ = t;
}
if (n==EulParOdd)
{
ea.fX = -ea.fX; ea.fY = -ea.fY; ea.fZ = -ea.fZ;
}
ti = ea.fX; tj = ea.fY; th = ea.fZ;
ci = cos(ti); cj = cos(tj); ch = cos(th);
si = sin(ti); sj = sin(tj); sh = sin(th);
cc = ci*ch; cs = ci*sh; sc = si*ch; ss = si*sh;
if (s==EulRepYes)
{
mat->fMap[i][i] = static_cast<float>(cj);
mat->fMap[i][j] = static_cast<float>(sj*si);
mat->fMap[i][k] = static_cast<float>(sj*ci);
mat->fMap[j][i] = static_cast<float>(sj*sh);
mat->fMap[j][j] = static_cast<float>(-cj*ss+cc);
mat->fMap[j][k] = static_cast<float>(-cj*cs-sc);
mat->fMap[k][i] = static_cast<float>(-sj*ch);
mat->fMap[k][j] = static_cast<float>(cj*sc+cs);
mat->fMap[k][k] = static_cast<float>(cj*cc-ss);
}
else
{
mat->fMap[i][i] = static_cast<float>(cj*ch);
mat->fMap[i][j] = static_cast<float>(sj*sc-cs);
mat->fMap[i][k] = static_cast<float>(sj*cc+ss);
mat->fMap[j][i] = static_cast<float>(cj*sh);
mat->fMap[j][j] = static_cast<float>(sj*ss+cc);
mat->fMap[j][k] = static_cast<float>(sj*cs-sc);
mat->fMap[k][i] = static_cast<float>(-sj);
mat->fMap[k][j] = static_cast<float>(cj*si);
mat->fMap[k][k] = static_cast<float>(cj*ci);
}
mat->fMap[W][X]=mat->fMap[W][Y]=mat->fMap[W][Z]=mat->fMap[X][W]=mat->fMap[Y][W]=mat->fMap[Z][W]=0.0;
mat->fMap[W][W]=1.0;
}
//
// Convert matrix to Euler angles (in radians)
//
void hsEuler::SetFromMatrix44(const hsMatrix44* mat, UInt32 order)
{
int i,j,k,h,n,s,f;
EulGetOrd(order,i,j,k,h,n,s,f);
if (s==EulRepYes)
{
double sy = sqrt(mat->fMap[i][j]*mat->fMap[i][j] + mat->fMap[i][k]*mat->fMap[i][k]);
if (sy > 16*FLT_EPSILON)
{
fX = static_cast<float>(atan2(mat->fMap[i][j], mat->fMap[i][k]));
fY = static_cast<float>(atan2(sy, (double)mat->fMap[i][i]));
fZ = static_cast<float>(atan2(mat->fMap[j][i], -mat->fMap[k][i]));
} else
{
fX = static_cast<float>(atan2(-mat->fMap[j][k], mat->fMap[j][j]));
fY = static_cast<float>(atan2(sy, (double)mat->fMap[i][i]));
fZ = 0;
}
}
else
{
double cy = sqrt(mat->fMap[i][i]*mat->fMap[i][i] + mat->fMap[j][i]*mat->fMap[j][i]);
if (cy > 16*FLT_EPSILON)
{
fX = static_cast<float>(atan2(mat->fMap[k][j], mat->fMap[k][k]));
fY = static_cast<float>(atan2((double)(-mat->fMap[k][i]), cy));
fZ = static_cast<float>(atan2(mat->fMap[j][i], mat->fMap[i][i]));
}
else
{
fX = static_cast<float>(atan2(-mat->fMap[j][k], mat->fMap[j][j]));
fY = static_cast<float>(atan2((double)(-mat->fMap[k][i]), cy));
fZ = 0;
}
}
if (n==EulParOdd)
{
fX = -fX; fY = - fY; fZ = -fZ;
}
if (f==EulFrmR)
{
hsScalar t = fX; fX = fZ; fZ = t;
}
fOrder = order;
}
//
// Convert quaternion to Euler angles (in radians)
//
void hsEuler::SetFromQuat(const hsQuat* q, UInt32 order)
{
hsMatrix44 mat;
double Nq = q->fX*q->fX+q->fY*q->fY+q->fZ*q->fZ+q->fW*q->fW;
double s = (Nq > 0.0) ? (2.0 / Nq) : 0.0;
double xs = q->fX*s, ys = q->fY*s, zs = q->fZ*s;
double wx = q->fW*xs, wy = q->fW*ys, wz = q->fW*zs;
double xx = q->fX*xs, xy = q->fX*ys, xz = q->fX*zs;
double yy = q->fY*ys, yz = q->fY*zs, zz = q->fZ*zs;
mat.fMap[X][X] = static_cast<float>(1.0 - (yy + zz));
mat.fMap[X][Y] = static_cast<float>(xy - wz);
mat.fMap[X][Z] = static_cast<float>(xz + wy);
mat.fMap[Y][X] = static_cast<float>(xy + wz);
mat.fMap[Y][Y] = static_cast<float>(1.0 - (xx + zz));
mat.fMap[Y][Z] = static_cast<float>(yz - wx);
mat.fMap[Z][X] = static_cast<float>(xz - wy);
mat.fMap[Z][Y] = static_cast<float>(yz + wx);
mat.fMap[Z][Z] = static_cast<float>(1.0 - (xx + yy));
mat.fMap[W][X] = mat.fMap[W][Y] = mat.fMap[W][Z] =
mat.fMap[X][W] = mat.fMap[Y][W] = mat.fMap[Z][W] = 0.0;
mat.fMap[W][W] = 1.0;
SetFromMatrix44(&mat, order);
}

116
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/hsEuler.h Normal file
View File

@ -0,0 +1,116 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
#ifndef HS_EULER_inc
#define HS_EULER_inc
#include "hsGeometry3.h"
//
//////////////////////////////////////////////////////////////////////////
// EULER STUFF
// See Gems IV, Ken Shoemake
//////////////////////////////////////////////////////////////////////////
//
/*** Order type constants, constructors, extractors ***/
/* There are 24 possible conventions, designated by: */
/* o EulAxI = axis used initially */
/* o EulPar = parity of axis permutation */
/* o EulRep = repetition of initial axis as last */
/* o EulFrm = frame from which axes are taken */
/* Axes I,J,K will be a permutation of X,Y,Z. */
/* Axis H will be either I or K, depending on EulRep. */
/* Frame S takes axes from initial static frame. */
/* If ord = (AxI=X, Par=Even, Rep=No, Frm=S), then */
/* {a,b,c,ord} means Rz(c)Ry(b)Rx(a), where Rz(c)v */
/* rotates v around Z by c radians. */
#define EulFrmS 0
#define EulFrmR 1
#define EulFrm(ord) ((unsigned)(ord)&1)
#define EulRepNo 0
#define EulRepYes 1
#define EulRep(ord) (((unsigned)(ord)>>1)&1)
#define EulParEven 0
#define EulParOdd 1
#define EulPar(ord) (((unsigned)(ord)>>2)&1)
#define EulSafe "\000\001\002\000"
#define EulNext "\001\002\000\001"
#define EulAxI(ord) ((int)(EulSafe[(((unsigned)(ord)>>3)&3)]))
#define EulAxJ(ord) ((int)(EulNext[EulAxI(ord)+(EulPar(ord)==EulParOdd)]))
#define EulAxK(ord) ((int)(EulNext[EulAxI(ord)+(EulPar(ord)!=EulParOdd)]))
#define EulAxH(ord) ((EulRep(ord)==EulRepNo)?EulAxK(ord):EulAxI(ord))
/* EulGetOrd unpacks all useful information about order simultaneously. */
#define EulGetOrd(ord,i,j,k,h,n,s,f) {unsigned o=ord;f=o&1;o>>=1;s=o&1;o>>=1;\
n=o&1;o>>=1;i=EulSafe[o&3];j=EulNext[i+n];k=EulNext[i+1-n];h=s?k:i;}
/* EulOrd creates an order value between 0 and 23 from 4-tuple choices. */
#define EulOrd(i,p,r,f) (((((((i)<<1)+(p))<<1)+(r))<<1)+(f))
/* Static axes */
#define EulOrdXYZs EulOrd(X,EulParEven,EulRepNo,EulFrmS)
#define EulOrdXYXs EulOrd(X,EulParEven,EulRepYes,EulFrmS)
#define EulOrdXZYs EulOrd(X,EulParOdd,EulRepNo,EulFrmS)
#define EulOrdXZXs EulOrd(X,EulParOdd,EulRepYes,EulFrmS)
#define EulOrdYZXs EulOrd(Y,EulParEven,EulRepNo,EulFrmS)
#define EulOrdYZYs EulOrd(Y,EulParEven,EulRepYes,EulFrmS)
#define EulOrdYXZs EulOrd(Y,EulParOdd,EulRepNo,EulFrmS)
#define EulOrdYXYs EulOrd(Y,EulParOdd,EulRepYes,EulFrmS)
#define EulOrdZXYs EulOrd(Z,EulParEven,EulRepNo,EulFrmS)
#define EulOrdZXZs EulOrd(Z,EulParEven,EulRepYes,EulFrmS)
#define EulOrdZYXs EulOrd(Z,EulParOdd,EulRepNo,EulFrmS)
#define EulOrdZYZs EulOrd(Z,EulParOdd,EulRepYes,EulFrmS)
/* Rotating axes */
#define EulOrdZYXr EulOrd(X,EulParEven,EulRepNo,EulFrmR)
#define EulOrdXYXr EulOrd(X,EulParEven,EulRepYes,EulFrmR)
#define EulOrdYZXr EulOrd(X,EulParOdd,EulRepNo,EulFrmR)
#define EulOrdXZXr EulOrd(X,EulParOdd,EulRepYes,EulFrmR)
#define EulOrdXZYr EulOrd(Y,EulParEven,EulRepNo,EulFrmR)
#define EulOrdYZYr EulOrd(Y,EulParEven,EulRepYes,EulFrmR)
#define EulOrdZXYr EulOrd(Y,EulParOdd,EulRepNo,EulFrmR)
#define EulOrdYXYr EulOrd(Y,EulParOdd,EulRepYes,EulFrmR)
#define EulOrdYXZr EulOrd(Z,EulParEven,EulRepNo,EulFrmR)
#define EulOrdZXZr EulOrd(Z,EulParEven,EulRepYes,EulFrmR)
#define EulOrdXYZr EulOrd(Z,EulParOdd,EulRepNo,EulFrmR)
#define EulOrdZYZr EulOrd(Z,EulParOdd,EulRepYes,EulFrmR)
struct hsMatrix44;
class hsQuat;
class hsEuler
{
public:
hsScalar fX,fY,fZ;
UInt32 fOrder;
hsEuler(hsScalar ai, hsScalar aj, hsScalar ah, UInt32 order) : fX(ai),fY(aj),fZ(ah),fOrder(order) {}
// getters, converters
void GetQuat(hsQuat* res );
void GetMatrix44(hsMatrix44* M);
// setters, converters
void SetFromMatrix44(const hsMatrix44* M, UInt32 order);
void SetFromQuat(const hsQuat* q, UInt32 order);
};
#endif // HS_EULER_inc

536
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/mat_decomp.cpp Normal file
View File

@ -0,0 +1,536 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
/**** Decompose.c ****/
/* Ken Shoemake, 1993 */
//
// Gems IV. Polar Decomp
//
#include <math.h>
#include "mat_decomp.h"
#ifdef __MWERKS__
//#pragma optimization_level 0
#endif
/******* Matrix Preliminaries *******/
/** Fill out 3x3 matrix to 4x4 **/
#define mat_pad(A) (A[W][X]=A[X][W]=A[W][Y]=A[Y][W]=A[W][Z]=A[Z][W]=0,A[W][W]=1)
/** Copy nxn matrix A to C using "gets" for assignment **/
#define mat_copy(C,gets,A,n) {int i,j; for(i=0;i<n;i++) for(j=0;j<n;j++)\
C[i][j] gets (A[i][j]);}
/** Copy transpose of nxn matrix A to C using "gets" for assignment **/
#define mat_tpose(AT,gets,A,n) {int i,j; for(i=0;i<n;i++) for(j=0;j<n;j++)\
AT[i][j] gets (A[j][i]);}
/** Assign nxn matrix C the element-wise combination of A and B using "op" **/
#define mat_binop(C,gets,A,op,B,n) {int i,j; for(i=0;i<n;i++) for(j=0;j<n;j++)\
C[i][j] gets (A[i][j]) op (B[i][j]);}
/** Multiply the upper left 3x3 parts of A and B to get AB **/
void mat_mult(const HMatrix A, const HMatrix B, HMatrix AB)
{
int i, j;
for (i=0; i<3; i++) for (j=0; j<3; j++)
AB[i][j] = A[i][0]*B[0][j] + A[i][1]*B[1][j] + A[i][2]*B[2][j];
}
/** Return dot product of length 3 vectors va and vb **/
float vdot(float *va, float *vb)
{
return (va[0]*vb[0] + va[1]*vb[1] + va[2]*vb[2]);
}
/** Set v to cross product of length 3 vectors va and vb **/
void vcross(float *va, float *vb, float *v)
{
v[0] = va[1]*vb[2] - va[2]*vb[1];
v[1] = va[2]*vb[0] - va[0]*vb[2];
v[2] = va[0]*vb[1] - va[1]*vb[0];
}
/** Set MadjT to transpose of inverse of M times determinant of M **/
void adjoint_transpose(HMatrix M, HMatrix MadjT)
{
vcross(M[1], M[2], MadjT[0]);
vcross(M[2], M[0], MadjT[1]);
vcross(M[0], M[1], MadjT[2]);
}
/******* Quaternion Preliminaries *******/
/* Construct a (possibly non-unit) quaternion from real components. */
gemQuat Qt_(float x, float y, float z, float w)
{
gemQuat qq;
qq.x = x; qq.y = y; qq.z = z; qq.w = w;
return (qq);
}
/* Return conjugate of quaternion. */
gemQuat Qt_Conj(gemQuat q)
{
gemQuat qq;
qq.x = -q.x; qq.y = -q.y; qq.z = -q.z; qq.w = q.w;
return (qq);
}
/* Return quaternion product qL * qR. Note: order is important!
* To combine rotations, use the product Mul(qSecond, qFirst),
* which gives the effect of rotating by qFirst then qSecond. */
gemQuat Qt_Mul(gemQuat qL, gemQuat qR)
{
gemQuat qq;
qq.w = qL.w*qR.w - qL.x*qR.x - qL.y*qR.y - qL.z*qR.z;
qq.x = qL.w*qR.x + qL.x*qR.w + qL.y*qR.z - qL.z*qR.y;
qq.y = qL.w*qR.y + qL.y*qR.w + qL.z*qR.x - qL.x*qR.z;
qq.z = qL.w*qR.z + qL.z*qR.w + qL.x*qR.y - qL.y*qR.x;
return (qq);
}
/* Return product of quaternion q by scalar w. */
gemQuat Qt_Scale(gemQuat q, float w)
{
gemQuat qq;
qq.w = q.w*w; qq.x = q.x*w; qq.y = q.y*w; qq.z = q.z*w;
return (qq);
}
/* Construct a unit quaternion from rotation matrix. Assumes matrix is
* used to multiply column vector on the left: vnew = mat vold. Works
* correctly for right-handed coordinate system and right-handed rotations.
* Translation and perspective components ignored. */
gemQuat Qt_FromMatrix(HMatrix mat)
{
/* This algorithm avoids near-zero divides by looking for a large component
* - first w, then x, y, or z. When the trace is greater than zero,
* |w| is greater than 1/2, which is as small as a largest component can be.
* Otherwise, the largest diagonal entry corresponds to the largest of |x|,
* |y|, or |z|, one of which must be larger than |w|, and at least 1/2. */
gemQuat qu;
register double tr, s;
tr = mat[X][X] + mat[Y][Y]+ mat[Z][Z];
if (tr >= 0.0) {
s = sqrt(tr + mat[W][W]);
qu.w = static_cast<float>(s*0.5);
s = 0.5 / s;
qu.x = static_cast<float>((mat[Z][Y] - mat[Y][Z]) * s);
qu.y = static_cast<float>((mat[X][Z] - mat[Z][X]) * s);
qu.z = static_cast<float>((mat[Y][X] - mat[X][Y]) * s);
} else {
int h = X;
if (mat[Y][Y] > mat[X][X]) h = Y;
if (mat[Z][Z] > mat[h][h]) h = Z;
switch (h) {
#define caseMacro(i,j,k,I,J,K) \
case I:\
s = sqrt( (mat[I][I] - (mat[J][J]+mat[K][K])) + mat[W][W] );\
qu.i = static_cast<float>(s*0.5);\
s = 0.5 / s;\
qu.j = static_cast<float>((mat[I][J] + mat[J][I]) * s);\
qu.k = static_cast<float>((mat[K][I] + mat[I][K]) * s);\
qu.w = static_cast<float>((mat[K][J] - mat[J][K]) * s);\
break
caseMacro(x,y,z,X,Y,Z);
caseMacro(y,z,x,Y,Z,X);
caseMacro(z,x,y,Z,X,Y);
}
}
if (mat[W][W] != 1.0) qu = Qt_Scale(qu, static_cast<float>(1/sqrt(mat[W][W])));
return (qu);
}
/******* Decomp Auxiliaries *******/
static HMatrix mat_id = {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
/** Compute either the 1 or infinity norm of M, depending on tpose **/
float mat_norm(HMatrix M, int tpose)
{
int i;
float sum, max;
max = 0.0;
for (i=0; i<3; i++) {
if (tpose) sum = static_cast<float>(fabs(M[0][i])+fabs(M[1][i])+fabs(M[2][i]));
else sum = static_cast<float>(fabs(M[i][0])+fabs(M[i][1])+fabs(M[i][2]));
if (max<sum) max = sum;
}
return max;
}
float norm_inf(HMatrix M) {return mat_norm(M, 0);}
float norm_one(HMatrix M) {return mat_norm(M, 1);}
/** Return index of column of M containing maximum abs entry, or -1 if M=0 **/
int find_max_col(HMatrix M)
{
float abs, max;
int i, j, col;
max = 0.0; col = -1;
for (i=0; i<3; i++) for (j=0; j<3; j++) {
abs = M[i][j]; if (abs<0.0) abs = -abs;
if (abs>max) {max = abs; col = j;}
}
return col;
}
/** Setup u for Household reflection to zero all v components but first **/
void make_reflector(float *v, float *u)
{
float s = static_cast<float>(sqrt(vdot(v, v)));
u[0] = v[0]; u[1] = v[1];
u[2] = v[2] + ((v[2]<0.0) ? -s : s);
s = static_cast<float>(sqrt(2.0/vdot(u, u)));
u[0] = u[0]*s; u[1] = u[1]*s; u[2] = u[2]*s;
}
/** Apply Householder reflection represented by u to column vectors of M **/
void reflect_cols(HMatrix M, float *u)
{
int i, j;
for (i=0; i<3; i++) {
float s = u[0]*M[0][i] + u[1]*M[1][i] + u[2]*M[2][i];
for (j=0; j<3; j++) M[j][i] -= u[j]*s;
}
}
/** Apply Householder reflection represented by u to row vectors of M **/
void reflect_rows(HMatrix M, float *u)
{
int i, j;
for (i=0; i<3; i++) {
float s = vdot(u, M[i]);
for (j=0; j<3; j++) M[i][j] -= u[j]*s;
}
}
/** Find orthogonal factor Q of rank 1 (or less) M **/
void do_rank1(HMatrix M, HMatrix Q)
{
float v1[3], v2[3], s;
int col;
mat_copy(Q,=,mat_id,4);
/* If rank(M) is 1, we should find a non-zero column in M */
col = find_max_col(M);
if (col<0) return; /* Rank is 0 */
v1[0] = M[0][col]; v1[1] = M[1][col]; v1[2] = M[2][col];
make_reflector(v1, v1); reflect_cols(M, v1);
v2[0] = M[2][0]; v2[1] = M[2][1]; v2[2] = M[2][2];
make_reflector(v2, v2); reflect_rows(M, v2);
s = M[2][2];
if (s<0.0) Q[2][2] = -1.0;
reflect_cols(Q, v1); reflect_rows(Q, v2);
}
/** Find orthogonal factor Q of rank 2 (or less) M using adjoint transpose **/
void do_rank2(HMatrix M, HMatrix MadjT, HMatrix Q)
{
float v1[3], v2[3];
float w, x, y, z, c, s, d;
int col;
/* If rank(M) is 2, we should find a non-zero column in MadjT */
col = find_max_col(MadjT);
if (col<0) {do_rank1(M, Q); return;} /* Rank<2 */
v1[0] = MadjT[0][col]; v1[1] = MadjT[1][col]; v1[2] = MadjT[2][col];
make_reflector(v1, v1); reflect_cols(M, v1);
vcross(M[0], M[1], v2);
make_reflector(v2, v2); reflect_rows(M, v2);
w = M[0][0]; x = M[0][1]; y = M[1][0]; z = M[1][1];
if (w*z>x*y) {
c = z+w; s = y-x; d = static_cast<float>(sqrt(c*c+s*s)); c = c/d; s = s/d;
Q[0][0] = Q[1][1] = c; Q[0][1] = -(Q[1][0] = s);
} else {
c = z-w; s = y+x; d = static_cast<float>(sqrt(c*c+s*s)); c = c/d; s = s/d;
Q[0][0] = -(Q[1][1] = c); Q[0][1] = Q[1][0] = s;
}
Q[0][2] = Q[2][0] = Q[1][2] = Q[2][1] = 0.0; Q[2][2] = 1.0;
reflect_cols(Q, v1); reflect_rows(Q, v2);
}
/******* Polar Decomposition *******/
/* Polar Decomposition of 3x3 matrix in 4x4,
* M = QS. See Nicholas Higham and Robert S. Schreiber,
* Fast Polar Decomposition of An Arbitrary Matrix,
* Technical Report 88-942, October 1988,
* Department of Computer Science, Cornell University.
*/
float polar_decomp(const HMatrix M, HMatrix Q, HMatrix S)
{
#define TOL 1.0e-6
HMatrix Mk, MadjTk, Ek;
float det, M_one, M_inf, MadjT_one, MadjT_inf, E_one, gamma, g1, g2;
int i, j;
mat_tpose(Mk,=,M,3);
M_one = norm_one(Mk); M_inf = norm_inf(Mk);
do {
adjoint_transpose(Mk, MadjTk);
det = vdot(Mk[0], MadjTk[0]);
if (det==0.0) {do_rank2(Mk, MadjTk, Mk); break;}
MadjT_one = norm_one(MadjTk); MadjT_inf = norm_inf(MadjTk);
gamma = static_cast<float>(sqrt(sqrt((MadjT_one*MadjT_inf)/(M_one*M_inf))/fabs(det)));
g1 = gamma*0.5f;
g2 = 0.5f/(gamma*det);
mat_copy(Ek,=,Mk,3);
mat_binop(Mk,=,g1*Mk,+,g2*MadjTk,3);
mat_copy(Ek,-=,Mk,3);
E_one = norm_one(Ek);
M_one = norm_one(Mk); M_inf = norm_inf(Mk);
} while (E_one>(M_one*TOL));
mat_tpose(Q,=,Mk,3); mat_pad(Q);
mat_mult(Mk, M, S); mat_pad(S);
for (i=0; i<3; i++) for (j=i; j<3; j++)
S[i][j] = S[j][i] = 0.5f*(S[i][j]+S[j][i]);
return (det);
}
/******* Spectral Decomposition *******/
/* Compute the spectral decomposition of symmetric positive semi-definite S.
* Returns rotation in U and scale factors in result, so that if K is a diagonal
* matrix of the scale factors, then S = U K (U transpose). Uses Jacobi method.
* See Gene H. Golub and Charles F. Van Loan. Matrix Computations. Hopkins 1983.
*/
HVect spect_decomp(HMatrix S, HMatrix U)
{
HVect kv;
double Diag[3],OffD[3]; /* OffD is off-diag (by omitted index) */
double g,h,fabsh,fabsOffDi,t,theta,c,s,tau,ta,OffDq,a,b;
static char nxt[] = {Y,Z,X};
int sweep, i, j;
mat_copy(U,=,mat_id,4);
Diag[X] = S[X][X]; Diag[Y] = S[Y][Y]; Diag[Z] = S[Z][Z];
OffD[X] = S[Y][Z]; OffD[Y] = S[Z][X]; OffD[Z] = S[X][Y];
for (sweep=20; sweep>0; sweep--) {
float sm = static_cast<float>(fabs(OffD[X])+fabs(OffD[Y])+fabs(OffD[Z]));
if (sm==0.0) break;
for (i=Z; i>=X; i--) {
int p = nxt[i]; int q = nxt[p];
fabsOffDi = fabs(OffD[i]);
g = 100.0*fabsOffDi;
if (fabsOffDi>0.0) {
h = Diag[q] - Diag[p];
fabsh = fabs(h);
if (fabsh+g==fabsh) {
t = OffD[i]/h;
} else {
theta = 0.5*h/OffD[i];
t = 1.0/(fabs(theta)+sqrt(theta*theta+1.0));
if (theta<0.0) t = -t;
}
c = 1.0/sqrt(t*t+1.0); s = t*c;
tau = s/(c+1.0);
ta = t*OffD[i]; OffD[i] = 0.0;
Diag[p] -= ta; Diag[q] += ta;
OffDq = OffD[q];
OffD[q] -= s*(OffD[p] + tau*OffD[q]);
OffD[p] += s*(OffDq - tau*OffD[p]);
for (j=Z; j>=X; j--) {
a = U[j][p]; b = U[j][q];
U[j][p] -= static_cast<float>(s*(b + tau*a));
U[j][q] += static_cast<float>(s*(a - tau*b));
}
}
}
}
kv.x = static_cast<float>(Diag[X]);
kv.y = static_cast<float>(Diag[Y]);
kv.z = static_cast<float>(Diag[Z]);
kv.w = 1.0f;
return (kv);
}
/******* Spectral Axis Adjustment *******/
/* Given a unit quaternion, q, and a scale vector, k, find a unit quaternion, p,
* which permutes the axes and turns freely in the plane of duplicate scale
* factors, such that q p has the largest possible w component, i.e. the
* smallest possible angle. Permutes k's components to go with q p instead of q.
* See Ken Shoemake and Tom Duff. Matrix Animation and Polar Decomposition.
* Proceedings of Graphics Interface 1992. Details on p. 262-263.
*/
gemQuat snuggle(gemQuat q, HVect *k)
{
#define SQRTHALF (0.7071067811865475244f)
#define sgn(n,v) ((n)?-(v):(v))
#define swap(a,i,j) {a[3]=a[i]; a[i]=a[j]; a[j]=a[3];}
#define cycle(a,p) if (p) {a[3]=a[0]; a[0]=a[1]; a[1]=a[2]; a[2]=a[3];}\
else {a[3]=a[2]; a[2]=a[1]; a[1]=a[0]; a[0]=a[3];}
gemQuat p;
float ka[4];
int i, turn = -1;
ka[X] = k->x; ka[Y] = k->y; ka[Z] = k->z;
if (ka[X]==ka[Y]) {if (ka[X]==ka[Z]) turn = W; else turn = Z;}
else {if (ka[X]==ka[Z]) turn = Y; else if (ka[Y]==ka[Z]) turn = X;}
if (turn>=0) {
gemQuat qtoz, qp;
unsigned neg[3], win;
double mag[3], t;
static gemQuat qxtoz = {0,SQRTHALF,0,SQRTHALF};
static gemQuat qytoz = {SQRTHALF,0,0,SQRTHALF};
static gemQuat qppmm = { 0.5, 0.5,-0.5,-0.5};
static gemQuat qpppp = { 0.5, 0.5, 0.5, 0.5};
static gemQuat qmpmm = {-0.5, 0.5,-0.5,-0.5};
static gemQuat qpppm = { 0.5, 0.5, 0.5,-0.5};
static gemQuat q0001 = { 0.0, 0.0, 0.0, 1.0};
static gemQuat q1000 = { 1.0, 0.0, 0.0, 0.0};
switch (turn) {
default: return (Qt_Conj(q));
case X: q = Qt_Mul(q, qtoz = qxtoz); swap(ka,X,Z) break;
case Y: q = Qt_Mul(q, qtoz = qytoz); swap(ka,Y,Z) break;
case Z: qtoz = q0001; break;
}
q = Qt_Conj(q);
mag[0] = (double)q.z*q.z+(double)q.w*q.w-0.5;
mag[1] = (double)q.x*q.z-(double)q.y*q.w;
mag[2] = (double)q.y*q.z+(double)q.x*q.w;
for (i=0; i<3; i++) if (neg[i] = (mag[i]<0.0)) mag[i] = -mag[i];
if (mag[0]>mag[1]) {if (mag[0]>mag[2]) win = 0; else win = 2;}
else {if (mag[1]>mag[2]) win = 1; else win = 2;}
switch (win) {
case 0: if (neg[0]) p = q1000; else p = q0001; break;
case 1: if (neg[1]) p = qppmm; else p = qpppp; cycle(ka,0) break;
case 2: if (neg[2]) p = qmpmm; else p = qpppm; cycle(ka,1) break;
}
qp = Qt_Mul(q, p);
t = sqrt(mag[win]+0.5);
p = Qt_Mul(p, Qt_(0.0,0.0,static_cast<float>(-qp.z/t),static_cast<float>(qp.w/t)));
p = Qt_Mul(qtoz, Qt_Conj(p));
} else {
float qa[4], pa[4];
unsigned lo, hi, neg[4], par = 0;
double all, big, two;
qa[0] = q.x; qa[1] = q.y; qa[2] = q.z; qa[3] = q.w;
for (i=0; i<4; i++) {
pa[i] = 0.0;
if (neg[i] = (qa[i]<0.0)) qa[i] = -qa[i];
par ^= neg[i];
}
/* Find two largest components, indices in hi and lo */
if (qa[0]>qa[1]) lo = 0; else lo = 1;
if (qa[2]>qa[3]) hi = 2; else hi = 3;
if (qa[lo]>qa[hi]) {
if (qa[lo^1]>qa[hi]) {hi = lo; lo ^= 1;}
else {hi ^= lo; lo ^= hi; hi ^= lo;}
} else {if (qa[hi^1]>qa[lo]) lo = hi^1;}
all = (qa[0]+qa[1]+qa[2]+qa[3])*0.5;
two = (qa[hi]+qa[lo])*SQRTHALF;
big = qa[hi];
if (all>two) {
if (all>big) {/*all*/
{int i; for (i=0; i<4; i++) pa[i] = static_cast<float>(sgn(neg[i], 0.5));}
cycle(ka,par)
} else {/*big*/ pa[hi] = static_cast<float>(sgn(neg[hi],1.0));}
} else {
if (two>big) {/*two*/
pa[hi] = static_cast<float>(sgn(neg[hi],SQRTHALF));
pa[lo] = static_cast<float>(sgn(neg[lo], SQRTHALF));
if (lo>hi) {hi ^= lo; lo ^= hi; hi ^= lo;}
if (hi==W) {hi = "\001\002\000"[lo]; lo = 3-hi-lo;}
swap(ka,hi,lo)
} else {/*big*/ pa[hi] = static_cast<float>(sgn(neg[hi],1.0));}
}
p.x = -pa[0]; p.y = -pa[1]; p.z = -pa[2]; p.w = pa[3];
}
k->x = ka[X]; k->y = ka[Y]; k->z = ka[Z];
return (p);
}
/******* Decompose Affine Matrix *******/
/* Decompose 4x4 affine matrix A as TFRUK(U transpose), where t contains the
* translation components, q contains the rotation R, u contains U, k contains
* scale factors, and f contains the sign of the determinant.
* Assumes A transforms column vectors in right-handed coordinates.
* See Ken Shoemake and Tom Duff. Matrix Animation and Polar Decomposition.
* Proceedings of Graphics Interface 1992.
*/
void decomp_affine(const HMatrix A, gemAffineParts *parts)
{
HMatrix Q, S, U;
gemQuat p;
float det;
parts->t = Qt_(A[X][W], A[Y][W], A[Z][W], 0);
det = polar_decomp(A, Q, S);
if (det<0.0) {
mat_copy(Q,=,-Q,3);
parts->f = -1;
} else parts->f = 1;
parts->q = Qt_FromMatrix(Q);
parts->k = spect_decomp(S, U);
parts->u = Qt_FromMatrix(U);
p = snuggle(parts->u, &parts->k);
parts->u = Qt_Mul(parts->u, p);
}
/******* Invert Affine Decomposition *******/
/* Compute inverse of affine decomposition.
*/
void invert_affine(gemAffineParts *parts, gemAffineParts *inverse)
{
gemQuat t, p;
inverse->f = parts->f;
inverse->q = Qt_Conj(parts->q);
inverse->u = Qt_Mul(parts->q, parts->u);
inverse->k.x = static_cast<float>((parts->k.x==0.0) ? 0.0 : 1.0/parts->k.x);
inverse->k.y = static_cast<float>((parts->k.y==0.0) ? 0.0 : 1.0/parts->k.y);
inverse->k.z = static_cast<float>((parts->k.z==0.0) ? 0.0 : 1.0/parts->k.z);
inverse->k.w = parts->k.w;
t = Qt_(-parts->t.x, -parts->t.y, -parts->t.z, 0);
t = Qt_Mul(Qt_Conj(inverse->u), Qt_Mul(t, inverse->u));
t = Qt_(inverse->k.x*t.x, inverse->k.y*t.y, inverse->k.z*t.z, 0);
p = Qt_Mul(inverse->q, inverse->u);
t = Qt_Mul(p, Qt_Mul(t, Qt_Conj(p)));
inverse->t = (inverse->f>0.0) ? t : Qt_(-t.x, -t.y, -t.z, 0);
}

57
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/mat_decomp.h Normal file
View File

@ -0,0 +1,57 @@
/*==LICENSE==*
CyanWorlds.com Engine - MMOG client, server and tools
Copyright (C) 2011 Cyan Worlds, Inc.
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 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
You can contact Cyan Worlds, Inc. by email legal@cyan.com
or by snail mail at:
Cyan Worlds, Inc.
14617 N Newport Hwy
Mead, WA 99021
*==LICENSE==*/
#ifndef MATDECOMP_inc
#define MATDECOMP_inc
//
// Types
//
typedef struct {float x, y, z, w;} gemQuat; /* Quaternion */
enum QuatPart {X, Y, Z, W};
typedef gemQuat HVect; /* Homogeneous 3D vector */
typedef float HMatrix[4][4]; /* Right-handed, for column vectors */
typedef struct {
HVect t; /* Translation components */
gemQuat q; /* Essential rotation */
gemQuat u; /* Stretch rotation */
HVect k; /* Stretch factors */
float f; /* Sign of determinant */
} gemAffineParts;
//
// Funcs
//
float polar_decomp(const HMatrix M, HMatrix Q, HMatrix S);
HVect spect_decomp(const HMatrix S, HMatrix U);
gemQuat snuggle(gemQuat q, HVect *k);
void decomp_affine(const HMatrix A, gemAffineParts *parts);
void invert_affine(gemAffineParts *parts, gemAffineParts *inverse);
#endif

1
MOULOpenSourceClientPlugin/Plasma20/Sources/Plasma/PubUtilLib/plTransform/notes.txt Normal file
View File

@ -0,0 +1 @@
-Didn't include hsGTransform.cpp,h because it references hsSceneObject, and plTMController.