OU/CWE-ou-minkata
2
3
Fork 0
mirror of https://foundry.openuru.org/gitblit/r/CWE-ou-minkata.git synced 2025-07-14 02:27:40 -04:00
Code Releases Activity

CWE Directory Reorganization

Browse Source 
Rearrange directory structure of CWE to be loosely equivalent to
the H'uru Plasma repository.

Part 1: Movement of directories and files.
This commit is contained in:
rarified
2021-05-15 12:49:46 -06:00
parent c3f4a640a3
commit 96903e8dca
4002 changed files with 159 additions and 644 deletions
Download Patch File Download Diff File Expand all files Collapse all files

442
Sources/Plasma/PubUtilLib/plTransform/hsAffineParts.cpp Normal file
View File

@ -0,0 +1,442 @@
/*==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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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);
}

89
Sources/Plasma/PubUtilLib/plTransform/hsAffineParts.h Normal file
View File

@ -0,0 +1,89 @@
/*==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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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

228
Sources/Plasma/PubUtilLib/plTransform/hsEuler.cpp Normal file
View File

@ -0,0 +1,228 @@
/*==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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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);
}

132
Sources/Plasma/PubUtilLib/plTransform/hsEuler.h Normal file
View File

@ -0,0 +1,132 @@
/*==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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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

552
Sources/Plasma/PubUtilLib/plTransform/mat_decomp.cpp Normal file
View File

@ -0,0 +1,552 @@
/*==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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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);
}

73
Sources/Plasma/PubUtilLib/plTransform/mat_decomp.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/>.
Additional permissions under GNU GPL version 3 section 7
If you modify this Program, or any covered work, by linking or
combining it with any of RAD Game Tools Bink SDK, Autodesk 3ds Max SDK,
NVIDIA PhysX SDK, Microsoft DirectX SDK, OpenSSL library, Independent
JPEG Group JPEG library, Microsoft Windows Media SDK, or Apple QuickTime SDK
(or a modified version of those libraries),
containing parts covered by the terms of the Bink SDK EULA, 3ds Max EULA,
PhysX SDK EULA, DirectX SDK EULA, OpenSSL and SSLeay licenses, IJG
JPEG Library README, Windows Media SDK EULA, or QuickTime SDK EULA, the
licensors of this Program grant you additional
permission to convey the resulting work. Corresponding Source for a
non-source form of such a combination shall include the source code for
the parts of OpenSSL and IJG JPEG Library used as well as that of the covered
work.
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
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.