From b4f17b968ad91f24c071ce1a1d1babf89fa9377c Mon Sep 17 00:00:00 2001 From: castano Date: Mon, 29 Dec 2008 11:27:13 +0000 Subject: [PATCH] Merge changes from internal branch. - Add frustum class and bezier evaluation functions. - Add component accessors to vector. - Add matrix constructors. - Fix errors in sparse solvers. - Better robust orthogonalization. - Fix montecarlo distribution. --- src/nvmath/Basis.cpp | 13 ++++- src/nvmath/Basis.h | 4 +- src/nvmath/Bezier.cpp | 97 ++++++++++++++++++++++++++++++++ src/nvmath/Bezier.h | 33 +++++++++++ src/nvmath/CMakeLists.txt | 4 +- src/nvmath/Frustum.cpp | 76 +++++++++++++++++++++++++ src/nvmath/Frustum.h | 34 +++++++++++ src/nvmath/Matrix.h | 19 ++++++- src/nvmath/Montecarlo.cpp | 31 ++-------- src/nvmath/Montecarlo.h | 31 ++++++++-- src/nvmath/Solver.cpp | 35 ++++++++---- src/nvmath/Sparse.cpp | 69 ++++++++++++++++------- src/nvmath/Sparse.h | 2 + src/nvmath/TypeSerialization.cpp | 5 ++ src/nvmath/TypeSerialization.h | 4 +- src/nvmath/Vector.h | 19 +++++++ 16 files changed, 408 insertions(+), 68 deletions(-) create mode 100644 src/nvmath/Bezier.cpp create mode 100644 src/nvmath/Bezier.h create mode 100644 src/nvmath/Frustum.cpp create mode 100644 src/nvmath/Frustum.h diff --git a/src/nvmath/Basis.cpp b/src/nvmath/Basis.cpp index 19e28ed..76660bb 100644 --- a/src/nvmath/Basis.cpp +++ b/src/nvmath/Basis.cpp @@ -27,7 +27,8 @@ void Basis::orthonormalize(float epsilon /*= NV_EPSILON*/) tangent -= normal * dot(normal, tangent); tangent = ::normalize(tangent, epsilon); - bitangent -= normal * dot(normal, bitangent) + tangent * dot(tangent, bitangent); + bitangent -= normal * dot(normal, bitangent); + bitangent -= tangent * dot(tangent, bitangent); bitangent = ::normalize(bitangent, epsilon); } @@ -91,7 +92,15 @@ void Basis::robustOrthonormalize(float epsilon /*= NV_EPSILON*/) tangent = nv::normalize(tangent); bitangent = nv::normalize(bitangent); - Vector3 bisector = nv::normalize(tangent + bitangent); + Vector3 bisector; + if (length(tangent + bitangent) < epsilon) + { + bisector = tangent; + } + else + { + bisector = nv::normalize(tangent + bitangent); + } Vector3 axis = cross(bisector, normal); nvDebugCheck(isNormalized(axis, epsilon)); diff --git a/src/nvmath/Basis.h b/src/nvmath/Basis.h index dca0442..0626437 100644 --- a/src/nvmath/Basis.h +++ b/src/nvmath/Basis.h @@ -12,8 +12,10 @@ namespace nv /// Basis class to compute tangent space basis, ortogonalizations and to /// transform vectors from one space to another. - struct Basis + class Basis { + public: + /// Create a null basis. Basis() : tangent(0, 0, 0), bitangent(0, 0, 0), normal(0, 0, 0) {} diff --git a/src/nvmath/Bezier.cpp b/src/nvmath/Bezier.cpp new file mode 100644 index 0000000..f108b2d --- /dev/null +++ b/src/nvmath/Bezier.cpp @@ -0,0 +1,97 @@ +// This code is in the public domain -- castanyo@yahoo.es + +#include "Bezier.h" + + +using namespace nv; + + +static void deCasteljau(float u, Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, Vector3 * p) +{ + Vector3 q0 = lerp(p0, p1, u); + Vector3 q1 = lerp(p1, p2, u); + Vector3 q2 = lerp(p2, p3, u); + Vector3 r0 = lerp(q0, q1, u); + Vector3 r1 = lerp(q1, q2, u); + + *p = lerp(r0, r1, u); +} + +static void deCasteljau(float u, Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, Vector3 * p, Vector3 * dp) +{ + Vector3 q0 = lerp(p0, p1, u); + Vector3 q1 = lerp(p1, p2, u); + Vector3 q2 = lerp(p2, p3, u); + Vector3 r0 = lerp(q0, q1, u); + Vector3 r1 = lerp(q1, q2, u); + + *dp = r0 - r1; + *p = lerp(r0, r1, u); +} + + +void nv::evaluateCubicBezierPatch(float u, float v, + Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, + Vector3::Arg p4, Vector3::Arg p5, Vector3::Arg p6, Vector3::Arg p7, + Vector3::Arg p8, Vector3::Arg p9, Vector3::Arg pA, Vector3::Arg pB, + Vector3::Arg pC, Vector3::Arg pD, Vector3::Arg pE, Vector3::Arg pF, + Vector3 * pos, Vector3 * du, Vector3 * dv) +{ +#if 0 + Vector2 L0(1-u,1-v); + Vector2 L1(u,v); + + Vector2 Q0 = L0 * L0; + Vector2 Q1 = 2 * L0 * L1; + Vector2 Q2 = L1 * L1; + + Vector2 B0 = L0 * L0 * L0; + Vector2 B1 = 3 * L0 * L0 * L1; + Vector2 B2 = 3 * L1 * L1 * L0; + Vector2 B3 = L1 * L1 * L1; + + *pos = + (B0.x() * p0 + B1.x() * p1 + B2.x() * p2 + B3.x() * p3) * B0.y() + + (B0.x() * p4 + B1.x() * p5 + B2.x() * p6 + B3.x() * p7) * B1.y() + + (B0.x() * p8 + B1.x() * p9 + B2.x() * pA + B3.x() * pB) * B2.y() + + (B0.x() * pC + B1.x() * pD + B2.x() * pE + B3.x() * pF) * B3.y(); + + *du = + ((p0-p1) * B0.y() + (p4-p5) * B1.y() + (p8-p9) * B2.y() + (pC-pD) * B3.y()) * Q0.x() + + ((p1-p2) * B0.y() + (p5-p6) * B1.y() + (p9-pA) * B2.y() + (pD-pE) * B3.y()) * Q1.x() + + ((p2-p3) * B0.y() + (p6-p7) * B1.y() + (pA-pB) * B2.y() + (pE-pF) * B3.y()) * Q2.x(); + + *dv = + ((p0-p4) * B0.x() + (p1-p5) * B1.x() + (p2-p6) * B2.x() + (p3-p7) * B3.x()) * Q0.y() + + ((p4-p8) * B0.x() + (p5-p9) * B1.x() + (p6-pA) * B2.x() + (p7-pB) * B3.x()) * Q1.y() + + ((p8-pC) * B0.x() + (p9-pD) * B1.x() + (pA-pE) * B2.x() + (pB-pF) * B3.x()) * Q2.y(); +#else + Vector3 t0, t1, t2, t3; + Vector3 q0, q1, q2, q3; + + deCasteljau(u, p0, p1, p2, p3, &q0, &t0); + deCasteljau(u, p4, p5, p6, p7, &q1, &t1); + deCasteljau(u, p8, p9, pA, pB, &q2, &t2); + deCasteljau(u, pC, pD, pE, pF, &q3, &t3); + + deCasteljau(v, q0, q1, q2, q3, pos, dv); + deCasteljau(v, t0, t1, t2, t3, du); +#endif +} + + +void nv::degreeElevate(Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, + Vector3 * q0, Vector3 * q1, Vector3 * q2, Vector3 * q3, Vector3 * q4) +{ +} + +void nv::evaluateQuarticBezierTriangle(float u, float v, + Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, Vector3::Arg p4, + Vector3::Arg p5, Vector3::Arg p6, Vector3::Arg p7, Vector3::Arg p8, + Vector3::Arg p9, Vector3::Arg pA, Vector3::Arg pB, + Vector3::Arg pC, Vector3::Arg pD, + Vector3::Arg pE, + Vector3 * pos, Vector3 * du, Vector3 * dv) +{ +} + diff --git a/src/nvmath/Bezier.h b/src/nvmath/Bezier.h new file mode 100644 index 0000000..2efe18a --- /dev/null +++ b/src/nvmath/Bezier.h @@ -0,0 +1,33 @@ +// This code is in the public domain -- castanyo@yahoo.es + +#ifndef NV_MATH_BEZIER_H +#define NV_MATH_BEZIER_H + +#include +#include + +namespace nv +{ + + void evaluateCubicBezierPatch(float u, float v, + Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, + Vector3::Arg p4, Vector3::Arg p5, Vector3::Arg p6, Vector3::Arg p7, + Vector3::Arg p8, Vector3::Arg p9, Vector3::Arg pA, Vector3::Arg pB, + Vector3::Arg pC, Vector3::Arg pD, Vector3::Arg pE, Vector3::Arg pF, + Vector3 * pos, Vector3 * du, Vector3 * dv); + + void degreeElevate(Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, + Vector3 * q0, Vector3 * q1, Vector3 * q2, Vector3 * q3, Vector3 * q4); + + void evaluateQuarticBezierTriangle(float u, float v, + Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector3::Arg p3, Vector3::Arg p4, + Vector3::Arg p5, Vector3::Arg p6, Vector3::Arg p7, Vector3::Arg p8, + Vector3::Arg p9, Vector3::Arg pA, Vector3::Arg pB, + Vector3::Arg pC, Vector3::Arg pD, + Vector3::Arg pE, + Vector3 * pos, Vector3 * du, Vector3 * dv); + + +} // nv namespace + +#endif // NV_MATH_BEZIER_H diff --git a/src/nvmath/CMakeLists.txt b/src/nvmath/CMakeLists.txt index 84b5a3e..bfe9bbf 100644 --- a/src/nvmath/CMakeLists.txt +++ b/src/nvmath/CMakeLists.txt @@ -8,6 +8,7 @@ SET(MATH_SRCS Plane.h Plane.cpp Box.h Color.h + Frustum.h Frustum.cpp Montecarlo.h Montecarlo.cpp Random.h Random.cpp SphericalHarmonic.h SphericalHarmonic.cpp @@ -19,7 +20,8 @@ SET(MATH_SRCS Solver.h Solver.cpp KahanSum.h Half.h Half.cpp - Fitting.h Fitting.cpp) + Fitting.h Fitting.cpp + Bezier.h Bezier.cpp) INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/nvmath/Frustum.cpp b/src/nvmath/Frustum.cpp new file mode 100644 index 0000000..60f6607 --- /dev/null +++ b/src/nvmath/Frustum.cpp @@ -0,0 +1,76 @@ +// This code is in the public domain -- castanyo@yahoo.es + +#include "Frustum.h" +#include "Box.h" +#include "Matrix.h" + +namespace nv +{ + Frustum::Frustum(float fovy, float aspect, float near, float far) + { + // Define a frustum looking along the -ve z axis. + planes[NEAR] = Plane(Vector3(0,0,1), near); + planes[FAR] = Plane(Vector3(0,0,-1), far); + + const Vector3 origin(zero); + const float tanfy = tan(fovy/2.0f); + planes[TOP] = Plane(Vector3( 1.0f, 0.0f, tanfy * aspect), origin); + planes[BOTTOM] = Plane(Vector3(-1.0f, 0.0f, tanfy * aspect), origin); + + planes[RIGHT] = Plane(Vector3(0.0f, 1.0f, tanfy), origin); + planes[LEFT] = Plane(Vector3(0.0f, -1.0f, tanfy), origin); + + for(int p = 0; p < Frustum::PLANE_COUNT; ++p) + planes[p] = normalize(planes[p]); + } + + static void getAllCorners(const Box& boxy, Vector3 c[]) + { + c[0] = boxy.minCorner(); + c[1] = boxy.maxCorner(); + c[2] = Vector3(c[0].x(), c[1].y(), c[0].z()); + c[3] = Vector3(c[0].x(), c[1].y(), c[1].z()); + c[4] = Vector3(c[0].x(), c[0].y(), c[1].z()); + c[5] = Vector3(c[1].x(), c[1].y(), c[0].z()); + c[6] = Vector3(c[1].x(), c[0].y(), c[0].z()); + c[7] = Vector3(c[1].x(), c[0].y(), c[1].z()); + } + + bool Frustum::intersects(const Box& boxy) const + { + const int N_CORNERS = 8; + Vector3 boxCorners[N_CORNERS]; + getAllCorners(boxy, boxCorners); + + // test all 8 corners against the 6 sides + // if all points are behind 1 specific plane, we are out + // if we are in with all points, then we are fully in + for(int p = 0; p < Frustum::PLANE_COUNT; ++p) + { + int iInCount = N_CORNERS; + + for(int i = 0; i < N_CORNERS; ++i) + { + // test point [i] against the planes + if (distance(planes[p], boxCorners[i]) > 0.0f) + --iInCount; + } + + // were all the points outside of plane p? + if (iInCount == 0) + return false; + } + + return true; + } + + Frustum transformFrustum(const Matrix& m, const Frustum& f) + { + Frustum result; + + for (int i=0; i!=Frustum::PLANE_COUNT; ++i) + result.planes[i] = transformPlane(m, f.planes[i]); + + return result; + } +} diff --git a/src/nvmath/Frustum.h b/src/nvmath/Frustum.h new file mode 100644 index 0000000..5a331ad --- /dev/null +++ b/src/nvmath/Frustum.h @@ -0,0 +1,34 @@ +// This code is in the public domain -- castanyo@yahoo.es + +#ifndef NV_MATH_FRUSTUM_H +#define NV_MATH_FRUSTUM_H + +#include + +namespace nv +{ + class Box; + + class NVMATH_CLASS Frustum + { + public: + // Construct a frustum from typical view parameters. fovy has the same meaning as for glut_perspective_reshaper. + // The resulting frustum is centred at the origin, pointing along the z-axis. + Frustum() {} + Frustum(float fovy, float aspect, float near, float far); + + // Unlike some intersection methods, we don't bother to distinguish intersects + // from contains. A true result could indicate either. + bool intersects(const Box&) const; + + private: + friend Frustum transformFrustum(const Matrix&, const Frustum&); + + enum PlaneIndices { NEAR, LEFT, RIGHT, TOP, BOTTOM, FAR, PLANE_COUNT }; + Plane planes[PLANE_COUNT]; + }; + + Frustum transformFrustum(const Matrix&, const Frustum&); +} + +#endif diff --git a/src/nvmath/Matrix.h b/src/nvmath/Matrix.h index 4064816..8a76896 100644 --- a/src/nvmath/Matrix.h +++ b/src/nvmath/Matrix.h @@ -24,6 +24,8 @@ public: Matrix(zero_t); Matrix(identity_t); Matrix(const Matrix & m); + Matrix(Vector4::Arg v0, Vector4::Arg v1, Vector4::Arg v2, Vector4::Arg v3); + Matrix(const scalar m[]); // m is assumed to contain 16 elements scalar data(uint idx) const; scalar & data(uint idx); @@ -75,6 +77,21 @@ inline Matrix::Matrix(const Matrix & m) } } +inline Matrix::Matrix(Vector4::Arg v0, Vector4::Arg v1, Vector4::Arg v2, Vector4::Arg v3) +{ + m_data[ 0] = v0.x(); m_data[ 1] = v0.y(); m_data[ 2] = v0.z(); m_data[ 3] = v0.w(); + m_data[ 4] = v1.x(); m_data[ 5] = v1.y(); m_data[ 6] = v1.z(); m_data[ 7] = v1.w(); + m_data[ 8] = v2.x(); m_data[ 9] = v2.y(); m_data[10] = v2.z(); m_data[11] = v2.w(); + m_data[12] = v3.x(); m_data[13] = v3.y(); m_data[14] = v3.z(); m_data[15] = v3.w(); +} + +inline Matrix::Matrix(const scalar m[]) +{ + for(int i = 0; i < 16; i++) { + m_data[i] = m[i]; + } +} + // Accessors inline scalar Matrix::data(uint idx) const @@ -332,7 +349,7 @@ inline Matrix transpose(Matrix::Arg m) Matrix r; for (int i = 0; i < 4; i++) { - for (int j = 0; j < 4; i++) + for (int j = 0; j < 4; j++) { r(i, j) = m(j, i); } diff --git a/src/nvmath/Montecarlo.cpp b/src/nvmath/Montecarlo.cpp index 4cd23a5..1fddd13 100644 --- a/src/nvmath/Montecarlo.cpp +++ b/src/nvmath/Montecarlo.cpp @@ -30,15 +30,8 @@ void SampleDistribution::redistributeRandom(const Distribution dist) { float x = m_rand.getFloat(); float y = m_rand.getFloat(); - - // Map uniform distribution in the square to the (hemi)sphere. - if( dist == Distribution_Uniform ) { - m_sampleArray[i].setUV(acosf(1 - 2 * x), 2 * PI * y); - } - else { - nvDebugCheck(dist == Distribution_Cosine); - m_sampleArray[i].setUV(acosf(sqrtf(x)), 2 * PI * y); - } + + setSample(i, dist, x, y); } } @@ -55,15 +48,8 @@ void SampleDistribution::redistributeStratified(const Distribution dist) for(uint u = 0; u < sqrtSampleCount; u++, i++) { float x = (u + m_rand.getFloat()) / float(sqrtSampleCount); float y = (v + m_rand.getFloat()) / float(sqrtSampleCount); - - // Map uniform distribution in the square to the (hemi)sphere. - if( dist == Distribution_Uniform ) { - m_sampleArray[i].setUV(acosf(1 - 2 * x), 2 * PI * y); - } - else { - nvDebugCheck(dist == Distribution_Cosine); - m_sampleArray[i].setUV(acosf(sqrtf(x)), 2 * PI * y); - } + + setSample(i, dist, x, y); } } } @@ -141,14 +127,7 @@ void SampleDistribution::redistributeNRook(const Distribution dist) float x = (i + m_rand.getFloat()) / sampleCount; float y = (cells[i] + m_rand.getFloat()) / sampleCount; - // Map uniform distribution in the square to the (hemi)sphere. - if( dist == Distribution_Uniform ) { - m_sampleArray[i].setUV(acosf(1 - 2 * x), 2 * PI * y); - } - else { - nvDebugCheck(dist == Distribution_Cosine); - m_sampleArray[i].setUV(acosf(sqrtf(x)), 2 * PI * y); - } + setSample(i, dist, x, y); } delete [] cells; diff --git a/src/nvmath/Montecarlo.h b/src/nvmath/Montecarlo.h index efda315..20c8abc 100644 --- a/src/nvmath/Montecarlo.h +++ b/src/nvmath/Montecarlo.h @@ -23,23 +23,26 @@ public: // Distribution functions. enum Distribution { - Distribution_Uniform, - Distribution_Cosine + Distribution_UniformSphere, + Distribution_UniformHemisphere, + Distribution_CosineHemisphere }; /// Constructor. - SampleDistribution(int num) + SampleDistribution(uint num) { m_sampleArray.resize(num); } + + uint count() const { return m_sampleArray.count(); } - void redistribute(Method method=Method_NRook, Distribution dist=Distribution_Cosine); + void redistribute(Method method=Method_NRook, Distribution dist=Distribution_CosineHemisphere); /// Get parametric coordinates of the sample. - Vector2 sample(int i) { return m_sampleArray[i].uv; } + Vector2 sample(int i) const { return m_sampleArray[i].uv; } /// Get sample direction. - Vector3 sampleDir(int i) { return m_sampleArray[i].dir; } + Vector3 sampleDir(int i) const { return m_sampleArray[i].dir; } /// Get number of samples. uint sampleCount() const { return m_sampleArray.count(); } @@ -70,6 +73,22 @@ private: Vector2 uv; Vector3 dir; }; + + inline void setSample(uint i, Distribution dist, float x, float y) + { + // Map uniform distribution in the square to desired domain. + if( dist == Distribution_UniformSphere ) { + m_sampleArray[i].setUV(acosf(1 - 2 * x), 2 * PI * y); + } + else if( dist == Distribution_UniformHemisphere ) { + m_sampleArray[i].setUV(acosf(x), 2 * PI * y); + } + else { + nvDebugCheck(dist == Distribution_CosineHemisphere); + m_sampleArray[i].setUV(acosf(sqrtf(x)), 2 * PI * y); + } + } + /// Random seed. MTRand m_rand; diff --git a/src/nvmath/Solver.cpp b/src/nvmath/Solver.cpp index 785376b..b8eabbf 100644 --- a/src/nvmath/Solver.cpp +++ b/src/nvmath/Solver.cpp @@ -44,7 +44,15 @@ namespace void apply(const FullVector & x, FullVector & y) const { - y *= x; + nvDebugCheck(x.dimension() == m_inverseDiagonal.dimension()); + nvDebugCheck(y.dimension() == m_inverseDiagonal.dimension()); + + // @@ Wrap vector component-wise product into a separate function. + const uint D = x.dimension(); + for (uint i = 0; i < D; i++) + { + y[i] = m_inverseDiagonal[i] * x[i]; + } } private: @@ -69,11 +77,16 @@ void nv::LeastSquaresSolver(const SparseMatrix & A, const FullVector & b, FullVe const uint D = A.width(); + SparseMatrix At(A.height(), A.width()); + transpose(A, At); + FullVector Atb(D); - mult(Transposed, A, b, Atb); + //mult(Transposed, A, b, Atb); + mult(At, b, Atb); SparseMatrix AtA(D); - mult(Transposed, A, NoTransposed, A, AtA); + //mult(Transposed, A, NoTransposed, A, AtA); + mult(At, A, AtA); SymmetricSolver(AtA, Atb, x, epsilon); } @@ -182,10 +195,10 @@ void nv::SymmetricSolver(const SparseMatrix & A, const FullVector & b, FullVecto nvDebugCheck(A.height() == b.dimension()); nvDebugCheck(b.dimension() == x.dimension()); -// JacobiPreconditioner jacobi(A, true); + JacobiPreconditioner jacobi(A, true); + ConjugateGradientSolver(jacobi, A, b, x, epsilon); -// ConjugateGradientSolver(jacobi, A, b, x, epsilon); - ConjugateGradientSolver(A, b, x, epsilon); +// ConjugateGradientSolver(A, b, x, epsilon); } @@ -280,9 +293,9 @@ void nv::SymmetricSolver(const SparseMatrix & A, const FullVector & b, FullVecto beta = delta_new / delta_old; - // p = r + beta·p - copy(r, p); - saxpy(beta, p, r); + // p = beta·p + r + scal(beta, p); + saxpy(1, r, p); } return i; @@ -358,8 +371,8 @@ void nv::SymmetricSolver(const SparseMatrix & A, const FullVector & b, FullVecto beta = delta_new / delta_old; // p = s + beta·p - copy(s, p); - saxpy(beta, p, s); + scal(beta, p); + saxpy(1, s, p); } return i; diff --git a/src/nvmath/Sparse.cpp b/src/nvmath/Sparse.cpp index 757093d..47e7c71 100644 --- a/src/nvmath/Sparse.cpp +++ b/src/nvmath/Sparse.cpp @@ -1,4 +1,4 @@ -// This code is in the public domain -- Ignacio Castao +// This code is in the public domain -- Ignacio Castaņo #include #include @@ -235,7 +235,7 @@ void FullMatrix::madRow(uint y, float alpha, FullVector & v) const const uint count = v.dimension(); for (uint i = 0; i < count; i++) { - v[i] += alpha * m_array[y * count + i]; + v[i] += m_array[y * count + i]; } } @@ -432,8 +432,11 @@ void SparseMatrix::setCoefficient(uint x, uint y, float f) } } - Coefficient c = { x, f }; - m_array[y].append( c ); + if (f != 0.0f) + { + Coefficient c = { x, f }; + m_array[y].append( c ); + } } void SparseMatrix::addCoefficient(uint x, uint y, float f) @@ -441,18 +444,21 @@ void SparseMatrix::addCoefficient(uint x, uint y, float f) nvDebugCheck( x < width() ); nvDebugCheck( y < height() ); - const uint count = m_array[y].count(); - for (uint i = 0; i < count; i++) + if (f != 0.0f) { - if (m_array[y][i].x == x) + const uint count = m_array[y].count(); + for (uint i = 0; i < count; i++) { - m_array[y][i].v += f; - return; + if (m_array[y][i].x == x) + { + m_array[y][i].v += f; + return; + } } - } - Coefficient c = { x, f }; - m_array[y].append( c ); + Coefficient c = { x, f }; + m_array[y].append( c ); + } } void SparseMatrix::mulCoefficient(uint x, uint y, float f) @@ -470,8 +476,11 @@ void SparseMatrix::mulCoefficient(uint x, uint y, float f) } } - Coefficient c = { x, f }; - m_array[y].append( c ); + if (f != 0.0f) + { + Coefficient c = { x, f }; + m_array[y].append( c ); + } } @@ -753,6 +762,32 @@ static float dotColumnColumn(int y, const SparseMatrix & A, int x, const SparseM } +void nv::transpose(const SparseMatrix & A, SparseMatrix & B) +{ + nvDebugCheck(A.width() == B.height()); + nvDebugCheck(B.width() == A.height()); + + const uint w = A.width(); + for (uint x = 0; x < w; x++) + { + B.clearRow(x); + } + + const uint h = A.height(); + for (uint y = 0; y < h; y++) + { + const Array & row = A.getRow(y); + + const uint count = row.count(); + for (uint i = 0; i < count; i++) + { + const SparseMatrix::Coefficient & c = row[i]; + nvDebugCheck(c.x < w); + + B.setCoefficient(y, c.x, c.v); + } + } +} // C = A * B void nv::mult(const SparseMatrix & A, const SparseMatrix & B, SparseMatrix & C) @@ -814,13 +849,9 @@ void nv::sgemm(float alpha, Transpose TA, const SparseMatrix & A, Transpose TB, c += alpha * dotRowRow(y, A, x, B); } - if (c != 0.0f) - { - C.setCoefficient(x, y, c); - } + C.setCoefficient(x, y, c); } } - } // C = At * A diff --git a/src/nvmath/Sparse.h b/src/nvmath/Sparse.h index 29439f8..54a2cc8 100644 --- a/src/nvmath/Sparse.h +++ b/src/nvmath/Sparse.h @@ -175,6 +175,8 @@ namespace nv }; + NVMATH_API void transpose(const SparseMatrix & A, SparseMatrix & B); + NVMATH_API void mult(const SparseMatrix & M, const FullVector & x, FullVector & y); NVMATH_API void mult(Transpose TM, const SparseMatrix & M, const FullVector & x, FullVector & y); diff --git a/src/nvmath/TypeSerialization.cpp b/src/nvmath/TypeSerialization.cpp index d70ddeb..b55608c 100644 --- a/src/nvmath/TypeSerialization.cpp +++ b/src/nvmath/TypeSerialization.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include @@ -80,3 +81,7 @@ Stream & nv::operator<< (Stream & s, Box & box) return s << box.m_mins << box.m_maxs; } +Stream & nv::operator<< (Stream & s, Plane & plane) +{ + return s << plane.asVector(); +} diff --git a/src/nvmath/TypeSerialization.h b/src/nvmath/TypeSerialization.h index 3c50a24..ccda460 100644 --- a/src/nvmath/TypeSerialization.h +++ b/src/nvmath/TypeSerialization.h @@ -15,8 +15,9 @@ namespace nv class Matrix; class Quaternion; - struct Basis; + class Basis; class Box; + class Plane; NVMATH_API Stream & operator<< (Stream & s, Vector2 & obj); NVMATH_API Stream & operator<< (Stream & s, Vector3 & obj); @@ -26,6 +27,7 @@ namespace nv NVMATH_API Stream & operator<< (Stream & s, Quaternion & obj); NVMATH_API Stream & operator<< (Stream & s, Basis & obj); NVMATH_API Stream & operator<< (Stream & s, Box & obj); + NVMATH_API Stream & operator<< (Stream & s, Plane & obj); } // nv namespace diff --git a/src/nvmath/Vector.h b/src/nvmath/Vector.h index d474090..329531a 100644 --- a/src/nvmath/Vector.h +++ b/src/nvmath/Vector.h @@ -27,6 +27,7 @@ public: Vector2(Vector2::Arg v); const Vector2 & operator=(Vector2::Arg v); + void setComponent(uint idx, scalar f); scalar x() const; scalar y() const; @@ -116,6 +117,7 @@ public: const Vector3 & xyz() const; scalar component(uint idx) const; + void setComponent(uint idx, scalar f); const scalar * ptr() const; @@ -162,6 +164,14 @@ inline scalar Vector2::component(uint idx) const return 0.0f; } +inline void Vector2::setComponent(uint idx, float f) +{ + nvDebugCheck(idx < 2); + if (idx == 0) m_x = f; + else if (idx == 1) m_y = f; +} + + inline const scalar * Vector2::ptr() const { return &m_x; @@ -362,6 +372,15 @@ inline scalar Vector4::component(uint idx) const return 0.0f; } +inline void Vector4::setComponent(uint idx, float f) +{ + nvDebugCheck(idx < 4); + if (idx == 0) m_x = f; + else if (idx == 1) m_y = f; + else if (idx == 2) m_z = f; + else if (idx == 3) m_w = f; +} + inline const scalar * Vector4::ptr() const { return &m_x;