diff --git a/libs/main/math/pimathmatrix.h b/libs/main/math/pimathmatrix.h index 9dd8ee7c..645c9e22 100644 --- a/libs/main/math/pimathmatrix.h +++ b/libs/main/math/pimathmatrix.h @@ -1227,7 +1227,7 @@ public: ++crow; smat.swapRows(i, crow); mtmp.swapRows(i, crow); - if (sv != 0) sv->swap(i, crow); + if (sv != 0) sv->swapElements(i, crow); } for (uint j = i + 1; j < _V2D::rows_; ++j) { mul = smat.element(i, j) / smat.element(i, i); diff --git a/libs/main/math/pimathvector.h b/libs/main/math/pimathvector.h index 09b8d07b..010b74e3 100644 --- a/libs/main/math/pimathvector.h +++ b/libs/main/math/pimathvector.h @@ -28,85 +28,151 @@ template class PIMathMatrixT; +#define PIMATHVECTOR_ZERO_CMP Type(1E-100) + /// Vector templated -#define PIMV_FOR(v, s) for (uint v = s; v < Size; ++v) +#define PIMV_FOR for (uint i = 0; i < Size; ++i) template class PIP_EXPORT PIMathVectorT { typedef PIMathVectorT _CVector; static_assert(std::is_arithmetic::value, "Type must be arithmetic"); - static_assert(Size > 0, "Size count must be > 0"); + static_assert(Size > 0, "Size must be > 0"); public: - PIMathVectorT() {resize();} - PIMathVectorT(const PIVector & val) {resize(); PIMV_FOR(i, 0) c[i] = val[i];} - PIMathVectorT(const _CVector & st, const _CVector & fn) {resize(); set(st, fn);} + PIMathVectorT(const Type & v = Type()) {PIMV_FOR c[i] = v;} + PIMathVectorT(const PIVector & val) { + assert(Size == val.size()); + PIMV_FOR c[i] = val[i]; + } + static _CVector fromTwoPoints(const _CVector & st, const _CVector & fn) { + _CVector tv; + PIMV_FOR tv[i] = fn[i] - st[i]; + return tv; + } uint size() const {return Size;} - _CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;} - _CVector & set(const _CVector & st, const _CVector & fn) {PIMV_FOR(i, 0) c[i] = fn[i] - st[i]; return *this;} - _CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;} - _CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;} - Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;} + _CVector & fill(const Type & v) {PIMV_FOR c[i] = v; return *this;} + _CVector & move(const Type & v) {PIMV_FOR c[i] += v; return *this;} + _CVector & move(const _CVector & v) {PIMV_FOR c[i] += v[i]; return *this;} + Type lengthSqr() const { + Type tv(0); + PIMV_FOR tv += c[i] * c[i]; + return tv; + } Type length() const {return sqrt(lengthSqr());} - Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;} - Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);} - Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);} + Type manhattanLength() const { + Type tv(0); + PIMV_FOR tv += piAbs(c[i]); + return tv; + } + Type angleCos(const _CVector & v) const { + Type tv = v.length() * length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + return dot(v) / tv; + } + Type angleSin(const _CVector & v) const { + Type tv = angleCos(v); + return sqrt(Type(1) - tv * tv); + } Type angleRad(const _CVector & v) const {return acos(angleCos(v));} - Type angleDeg(const _CVector & v) const {return toDeg(acos(angleCos(v)));} - Type angleElevation(const _CVector & v) const {_CVector z = v - *this; double c = z.angleCos(*this); return 90.0 - acos(c) * rad2deg;} - _CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));} - _CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; if (piAbs(tv) <= Type(1E-100)) {fill(Type(0)); return *this;} PIMV_FOR(i, 0) c[i] /= tv; return *this;} + Type angleDeg(const _CVector & v) const {return toDeg(angleRad(v));} + Type angleElevation(const _CVector & v) const {return 90.0 - angleDeg(v - *this);} + _CVector projection(const _CVector & v) { + Type tv = v.length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + return v * (dot(v) / tv); + } + _CVector & normalize() { + Type tv = length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + if (tv == Type(1)) return *this; + PIMV_FOR c[i] /= tv; + return *this; + } _CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;} - _CVector cross(const _CVector & v) {return (*this) * v;} - Type dot(const _CVector & v) const {return (*this) ^ v;} - bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;} + bool isNull() const {PIMV_FOR if (c[i] != Type(0)) return false; return true;} bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);} Type & at(uint index) {return c[index];} Type at(uint index) const {return c[index];} Type & operator [](uint index) {return c[index];} Type operator [](uint index) const {return c[index];} - _CVector & operator =(const _CVector & v) {memcpy(c, v.c, sizeof(Type) * Size); return *this;} - _CVector & operator =(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;} - bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;} + _CVector & operator =(const Type & v) {PIMV_FOR c[i] = v; return *this;} + bool operator ==(const _CVector & v) const {PIMV_FOR if (c[i] != v[i]) return false; return true;} bool operator !=(const _CVector & v) const {return !(*this == c);} - void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];} - void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];} - void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;} - void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];} - void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;} - void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];} - _CVector operator -() const {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;} - _CVector operator +(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;} - _CVector operator -(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;} - _CVector operator *(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;} - _CVector operator /(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;} - _CVector operator /(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v[i]; return tv;} - _CVector operator *(const _CVector & v) const {if (Size != 3) return _CVector(); _CVector tv; tv.fill(Type(1)); tv[0] = c[1]*v[2] - v[1]*c[2]; tv[1] = v[0]*c[2] - c[0]*v[2]; tv[2] = c[0]*v[1] - v[0]*c[1]; return tv;} - _CVector operator &(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v[i]; return tv;} - Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;} + void operator +=(const _CVector & v) {PIMV_FOR c[i] += v[i];} + void operator -=(const _CVector & v) {PIMV_FOR c[i] -= v[i];} + void operator *=(const Type & v) {PIMV_FOR c[i] *= v;} + void operator /=(const Type & v) { + assert(piAbs(v) > PIMATHVECTOR_ZERO_CMP); + PIMV_FOR c[i] /= v; + } + _CVector operator -() const { + _CVector tv; + PIMV_FOR tv[i] = -c[i]; + return tv; + } + _CVector operator +(const _CVector & v) const { + _CVector tv(*this); + PIMV_FOR tv[i] += v[i]; + return tv; + } + _CVector operator -(const _CVector & v) const { + _CVector tv(*this); + PIMV_FOR tv[i] -= v[i]; + return tv; + } + _CVector operator *(const Type & v) const { + _CVector tv(*this); + PIMV_FOR tv[i] *= v; + return tv; + } + _CVector operator /(const Type & v) const { + assert(piAbs(v) > PIMATHVECTOR_ZERO_CMP); + _CVector tv = _CVector(*this); + PIMV_FOR tv[i] /= v; + return tv; + } + + _CVector cross(const _CVector & v) const { + static_assert(Size == 3, "cross product avalible only for 3D vectors"); + _CVector tv; + tv[0] = c[1]*v[2] - v[1]*c[2]; + tv[1] = v[0]*c[2] - c[0]*v[2]; + tv[2] = c[0]*v[1] - v[0]*c[1]; + return tv; + } + Type dot(const _CVector & v) const { + Type tv(0); + PIMV_FOR tv += c[i] * v[i]; + return tv; + } PIMathMatrixT<1, Size, Type> transposed() const { PIMathMatrixT<1, Size, Type> ret; - PIMV_FOR(i, 0) ret[0][i] = c[i]; + PIMV_FOR ret[0][i] = c[i]; return ret; } Type distToLine(const _CVector & lp0, const _CVector & lp1) { - _CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this); - Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length(); - return f;} + _CVector a(lp0, lp1); + Type tv = a.length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + _CVector b(lp0, *this); + return piAbs(a[0]*b[1] - a[1]*b[0]) / tv; + } template /// vector {Size, Type} to vector {Size1, Type1} - PIMathVectorT turnTo() const {PIMathVectorT tv; uint sz = piMin(Size, Size1); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;} - - static _CVector filled(const Type & v) {_CVector vv; PIMV_FOR(i, 0) vv[i] = v; return vv;} + PIMathVectorT turnTo() const { + PIMathVectorT tv; + uint sz = piMin(Size, Size1); + for (uint i = 0; i < sz; ++i) tv[i] = c[i]; + return tv; + } private: - void resize(const Type & new_value = Type()) {for (uint i = 0; i < Size; ++i) c[i] = new_value;} - Type c[Size]; }; @@ -117,18 +183,8 @@ inline PIMathVectorT operator *(const Type & x, const PIMathVectorT< } template -inline PICout operator <<(PICout s, const PIMathVectorT & v) {s << "{"; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << "}"; return s;} -template -inline bool operator ||(const PIMathVectorT & f, const PIMathVectorT & s) {return (f * s).isNull();} -template -inline PIMathVectorT sqrt(const PIMathVectorT & v) {PIMathVectorT ret; PIMV_FOR(i, 0) {ret[i] = sqrt(v[i]);} return ret;} -template -inline PIMathVectorT sqr(const PIMathVectorT & v) {PIMathVectorT ret; PIMV_FOR(i, 0) {ret[i] = sqr(v[i]);} return ret;} +inline PICout operator <<(PICout s, const PIMathVectorT & v) {s << "{"; PIMV_FOR {s << v[i]; if (i < Size - 1) s << ", ";} s << "}"; return s;} -template -inline PIByteArray & operator <<(PIByteArray & s, const PIMathVectorT & v) {for (uint i = 0; i < Size; ++i) s << v[i]; return s;} -template -inline PIByteArray & operator >>(PIByteArray & s, PIMathVectorT & v) {for (uint i = 0; i < Size; ++i) s >> v[i]; return s;} template @@ -158,7 +214,7 @@ typedef PIMathVectorT<4u, double> PIMathVectorT4d; /// Vector -#define PIMV_FOR(v, s) for (uint v = s; v < c.size(); ++v) +#define PIMV_FOR for (uint i = 0; i < c.size(); ++i) template class PIP_EXPORT PIMathVector { @@ -166,63 +222,167 @@ class PIP_EXPORT PIMathVector { template friend PIByteArray & operator <<(PIByteArray & s, const PIMathVector & v); template friend PIByteArray & operator >>(PIByteArray & s, PIMathVector & v); public: - PIMathVector(const uint size = 0) {c.resize(size);} - PIMathVector(const PIVector & val) {c.resize(val.size()); PIMV_FOR(i, 0) c[i] = val[i];} - PIMathVector(const _CVector & st, const _CVector & fn) {c.resize(st.size()); PIMV_FOR(i, 0) c[i] = fn[i] - st[i];} + PIMathVector(const uint size = 0, const Type & new_value = Type()) {c.resize(size, new_value);} + PIMathVector(const PIVector & val) {c = val;} + + template + PIMathVector(const PIMathVectorT & val) {c.resize(Size); PIMV_FOR c[i] = val[i];} + + static PIMathVector fromTwoPoints(const _CVector & st, const _CVector & fn) { + assert(st.size() == fn.size()); + _CVector v(st.size()); + for (uint i = 0; i < v.size(); ++i) v.c[i] = fn[i] - st[i]; + } uint size() const {return c.size();} - _CVector & resize(uint size, const Type & new_value = Type()) {c.resize(size, new_value); return *this;} - _CVector resized(uint size, const Type & new_value = Type()) {_CVector tv = _CVector(*this); tv.resize(size, new_value); return tv;} - _CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;} - _CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;} - _CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;} - _CVector & swap(uint fe, uint se) {piSwap(c[fe], c[se]); return *this;} - Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;} + _CVector & resize(uint size, const Type & new_value = Type()) { + c.resize(size, new_value); + return *this; + } + _CVector resized(uint size, const Type & new_value = Type()) { + _CVector tv = _CVector(*this); + tv.resize(size, new_value); + return tv; + } + _CVector & fill(const Type & v) { + c.fill(v); + return *this; + } + _CVector & move(const Type & v) { + PIMV_FOR c[i] += v; + return *this; + } + _CVector & move(const _CVector & v) { + assert(c.size() == v.size()); + PIMV_FOR c[i] += v[i]; + return *this; + } + _CVector & swapElements(uint fe, uint se) { + piSwap(c[fe], c[se]); + return *this; + } + Type lengthSqr() const { + Type tv(0); + PIMV_FOR tv += c[i] * c[i]; + return tv; + } Type length() const {return sqrt(lengthSqr());} - Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;} - Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);} - Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);} + Type manhattanLength() const { + Type tv(0); + PIMV_FOR tv += piAbs(c[i]); + return tv; + } + Type angleCos(const _CVector & v) const { + assert(c.size() == v.size()); + Type tv = v.length() * length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + return dot(v) / tv; + } + Type angleSin(const _CVector & v) const { + assert(c.size() == v.size()); + Type tv = angleCos(v); + return sqrt(Type(1) - tv * tv); + } Type angleRad(const _CVector & v) const {return acos(angleCos(v));} - Type angleDeg(const _CVector & v) const {return toDeg(acos(angleCos(v)));} - _CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));} - _CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; if (piAbs(tv) <= Type(1E-100)) {fill(Type(0)); return *this;} PIMV_FOR(i, 0) c[i] /= tv; return *this;} - _CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;} - bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;} + Type angleDeg(const _CVector & v) const {return toDeg(angleRad(v));} + _CVector projection(const _CVector & v) { + assert(c.size() == v.size()); + Type tv = v.length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + return v * (dot(v) / tv); + } + _CVector & normalize() { + Type tv = length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + if (tv == Type(1)) return *this; + PIMV_FOR c[i] /= tv; + return *this; + } + _CVector normalized() { + _CVector tv(*this); + tv.normalize(); + return tv; + } + bool isNull() const { + PIMV_FOR if (c[i] != Type(0)) return false; + return true; + } bool isValid() const {return !c.isEmpty();} - - bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);} + bool isOrtho(const _CVector & v) const {return dot(v) == Type(0);} Type & at(uint index) {return c[index];} Type at(uint index) const {return c[index];} Type & operator [](uint index) {return c[index];} Type operator [](uint index) const {return c[index];} - _CVector & operator =(const _CVector & v) {c = v.c; return *this;} - _CVector & operator =(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;} - bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;} - bool operator !=(const _CVector & v) const {return !(*this == c);} - void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];} - void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];} - void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;} - void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];} - void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;} - void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];} - _CVector operator -() const {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;} - _CVector operator +(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;} - _CVector operator -(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;} - _CVector operator *(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;} - _CVector operator /(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;} - _CVector operator *(const _CVector & v) const {if (c.size() < 3) return _CVector(); _CVector tv; tv.fill(Type(1)); tv[0] = c[1]*v[2] - v[1]*c[2]; tv[1] = v[0]*c[2] - c[0]*v[2]; tv[2] = c[0]*v[1] - v[0]*c[1]; return tv;} - _CVector operator &(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v[i]; return tv;} - Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;} - - Type distToLine(const _CVector & lp0, const _CVector & lp1) { - _CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this); - Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length(); - return f; + bool operator ==(const _CVector & v) const {return c == v.c;} + bool operator !=(const _CVector & v) const {return c != v.c;} + void operator +=(const _CVector & v) { + assert(c.size() == v.size()); + PIMV_FOR c[i] += v[i]; + } + void operator -=(const _CVector & v) { + assert(c.size() == v.size()); + PIMV_FOR c[i] -= v[i]; + } + void operator *=(const Type & v) {PIMV_FOR c[i] *= v;} + void operator /=(const Type & v) { + assert(piAbs(v) > PIMATHVECTOR_ZERO_CMP); + PIMV_FOR c[i] /= v; + } + _CVector operator -() const { + _CVector tv(c.size()); + PIMV_FOR tv[i] = -c[i]; + return tv; + } + _CVector operator +(const _CVector & v) const { + assert(c.size() == v.size()); + _CVector tv(*this); + PIMV_FOR tv[i] += v[i]; + return tv; + } + _CVector operator -(const _CVector & v) const { + assert(c.size() == v.size()); + _CVector tv(*this); + PIMV_FOR tv[i] -= v[i]; + return tv; + } + _CVector operator *(const Type & v) const { + _CVector tv(*this); + PIMV_FOR tv[i] *= v; + return tv; + } + _CVector operator /(const Type & v) const { + assert(piAbs(v) > PIMATHVECTOR_ZERO_CMP); + _CVector tv(*this); + PIMV_FOR tv[i] /= v; + return tv; + } + _CVector cross(const _CVector & v) const { + assert(c.size() == 3); + assert(v.size() == 3); + _CVector tv(3); + tv[0] = c[1]*v[2] - v[1]*c[2]; + tv[1] = c[2]*v[0] - v[2]*c[0]; + tv[2] = c[0]*v[1] - v[0]*c[1]; + return tv; + } + Type dot(const _CVector & v) const { + assert(c.size() == v.size()); + Type tv(0); + PIMV_FOR tv += c[i] * v[i]; + return tv; + } + + Type distToLine(const _CVector & lp0, const _CVector & lp1) { + assert(c.size() == lp0.size()); + assert(c.size() == lp1.size()); + _CVector a = _CVector::fromTwoPoints(lp0, lp1); + Type tv = a.length(); + assert(piAbs(tv) > PIMATHVECTOR_ZERO_CMP); + _CVector b = _CVector::fromTwoPoints(lp0, *this); + return piAbs(a[0]*b[1] - a[1]*b[0]) / tv; } - template - PIMathVector turnTo(uint size) const {PIMathVector tv; uint sz = piMin(c.size(), size); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;} PIVector toVector() const {return c;} inline Type * data() {return c.data();} @@ -233,6 +393,11 @@ private: }; +template +inline PIMathVector operator *(const Type & x, const PIMathVector & v) { + return v * x; +} + #undef PIMV_FOR #ifdef PIP_STD_IOSTREAM diff --git a/libs/main/math/piquaternion.cpp b/libs/main/math/piquaternion.cpp index 5ea23fe9..fa95ad7c 100644 --- a/libs/main/math/piquaternion.cpp +++ b/libs/main/math/piquaternion.cpp @@ -199,8 +199,8 @@ PIQuaternion PIQuaternion::fromRotationMatrix(const PIMathMatrixT33d & m) { PIQuaternion operator*(const PIQuaternion & q0, const PIQuaternion & q1) { PIMathVectorT3d v0(q0.vector()), v1(q1.vector()); - double r0 = q0.q[0] * q1.q[0] - (v0^v1); - PIMathVectorT3d qv = v1*q0.q[0] + v0*q1.q[0] + v0*v1; + double r0 = q0.q[0] * q1.q[0] - v0.dot(v1); + PIMathVectorT3d qv = v1*q0.q[0] + v0*q1.q[0] + v0.cross(v1); PIQuaternion ret; ret.q[0] = r0; ret.q[1] = qv[0];