From 4c359f42a7826c009abd264e9f18cf4d7c10b736 Mon Sep 17 00:00:00 2001 From: Pavel Krajcevski Date: Mon, 8 Oct 2012 17:53:27 -0400 Subject: [PATCH] - Added a parameter to the PCA computation that returns the first and second eigenvalues of the covariance matrix associated with the cluster. - Compared results of testing the ratio of eigenvalues as a measurement of 'linearity' for the different shapes, and output statistics. - Added a #define that controls whether or not we do shape estimation using quantized AABB error or eigenvalue ratios. The former seems to be better. --- BPTCEncoder/src/BC7Compressor.cpp | 172 ++++++++++++++++++++++++++---- BPTCEncoder/src/RGBAEndpoints.cpp | 161 +++++++++++++++++++++------- BPTCEncoder/src/RGBAEndpoints.h | 25 ++++- 3 files changed, 297 insertions(+), 61 deletions(-) diff --git a/BPTCEncoder/src/BC7Compressor.cpp b/BPTCEncoder/src/BC7Compressor.cpp index 0a8de51..8a5ecdc 100755 --- a/BPTCEncoder/src/BC7Compressor.cpp +++ b/BPTCEncoder/src/BC7Compressor.cpp @@ -37,10 +37,16 @@ #define ALIGN_SSE __attribute__((aligned(16))) #endif +#define USE_PCA_FOR_SHAPE_ESTIMATION + enum EBlockStats { eBlockStat_Path, eBlockStat_Mode, + eBlockStat_SingleShapeEstimate, + eBlockStat_TwoShapeEstimate, + eBlockStat_ThreeShapeEstimate, + eBlockStat_ModeZeroEstimate, eBlockStat_ModeOneEstimate, eBlockStat_ModeTwoEstimate, @@ -66,6 +72,10 @@ static const char *kBlockStatString[kNumBlockStats] = { "BlockStat_Path", "BlockStat_Mode", + "BlockStat_SingleShapeEstimate", + "BlockStat_TwoShapeEstimate", + "BlockStat_ThreeShapeEstimate", + "BlockStat_ModeZeroEstimate", "BlockStat_ModeOneEstimate", "BlockStat_ModeTwoEstimate", @@ -707,7 +717,7 @@ double BC7CompressionMode::CompressCluster(const RGBACluster &cluster, RGBAVecto bestIndices[i] = 1; alphaIndices[i] = 1; } - + return bestErr; } @@ -972,7 +982,8 @@ double BC7CompressionMode::CompressCluster(const RGBACluster &cluster, RGBAVecto #if 1 RGBAVector avg = cluster.GetTotal() / float(cluster.GetNumPoints()); RGBADir axis; - ::GetPrincipalAxis(cluster.GetNumPoints(), cluster.GetPoints(), axis); + double eigOne; + ::GetPrincipalAxis(cluster.GetNumPoints(), cluster.GetPoints(), axis, eigOne, NULL); float mindp = FLT_MAX, maxdp = -FLT_MAX; for(int i = 0 ; i < cluster.GetNumPoints(); i++) { @@ -1677,55 +1688,91 @@ namespace BC7C assert(!(clusters[0].GetPointBitString() & clusters[2].GetPointBitString())); } - static double EstimateTwoClusterError(RGBACluster &c) { + static double EstimateTwoClusterError(RGBACluster &c, double (&estimates)[2]) { RGBAVector Min, Max, v; c.GetBoundingBox(Min, Max); v = Max - Min; if(v * v == 0) { - gModeEstimate[1] = gModeEstimate[3] = 0.0; + estimates[0] = estimates[1] = 0.0; return 0.0; } const float *w = BC7C::GetErrorMetric(); - const double err1 = 0.0001 + c.QuantizedError(Min, Max, 8, 0xFFFCFCFC, RGBAVector(w[0], w[1], w[2], w[3])); + const double err1 = c.QuantizedError(Min, Max, 8, 0xFFFCFCFC, RGBAVector(w[0], w[1], w[2], w[3])); if(err1 >= 0.0) - gModeEstimate[1] = err1; + estimates[0] = err1; else - gModeEstimate[1] = min(gModeEstimate[1], err1); + estimates[0] = min(estimates[0], err1); - const double err3 = 0.0001 + c.QuantizedError(Min, Max, 8, 0xFFFEFEFE, RGBAVector(w[0], w[1], w[2], w[3])); + const double err3 = c.QuantizedError(Min, Max, 8, 0xFFFEFEFE, RGBAVector(w[0], w[1], w[2], w[3])); if(err3 >= 0.0) - gModeEstimate[3] = err3; + estimates[1] = err3; else - gModeEstimate[3] = min(gModeEstimate[3], err3); + estimates[1] = min(estimates[1], err3); - return min(err1, err3); + double error = 0.0001; +#ifdef USE_PCA_FOR_SHAPE_ESTIMATION + double eigOne = c.GetPrincipalEigenvalue(); + double eigTwo = c.GetSecondEigenvalue(); + if(eigOne != 0.0) { + error += eigTwo / eigOne; + } + else { + error += 1.0; + } +#else + error += min(err1, err3); +#endif + return error; } - static double EstimateThreeClusterError(RGBACluster &c) { + static double EstimateThreeClusterError(RGBACluster &c, double (&estimates)[2]) { RGBAVector Min, Max, v; c.GetBoundingBox(Min, Max); v = Max - Min; if(v * v == 0) { - gModeEstimate[0] = gModeEstimate[2] = 0.0; + estimates[0] = estimates[1] = 0.0; return 0.0; } const float *w = BC7C::GetErrorMetric(); const double err0 = 0.0001 + c.QuantizedError(Min, Max, 4, 0xFFF0F0F0, RGBAVector(w[0], w[1], w[2], w[3])); if(err0 >= 0.0) - gModeEstimate[0] = err0; + estimates[0] = err0; else - gModeEstimate[0] = min(gModeEstimate[0], err0); + estimates[0] = min(estimates[0], err0); const double err2 = 0.0001 + c.QuantizedError(Min, Max, 4, 0xFFF8F8F8, RGBAVector(w[0], w[1], w[2], w[3])); if(err2 >= 0.0) - gModeEstimate[2] = err2; + estimates[1] = err2; else - gModeEstimate[2] = min(gModeEstimate[2], err2); + estimates[1] = min(estimates[1], err2); - return min(err0, err2); + double error = 0.0001; +#ifdef USE_PCA_FOR_SHAPE_ESTIMATION + double eigOne = c.GetPrincipalEigenvalue(); + double eigTwo = c.GetSecondEigenvalue(); + + // printf("EigOne: %08.3f\tEigTwo: %08.3f\n", eigOne, eigTwo); + if(eigOne != 0.0) { + error += eigTwo / eigOne; + } + else { + error += 1.0; + } +#else + error += min(err0, err2); +#endif + return error; + } + + static void UpdateErrorEstimate(uint32 mode, double est) { + assert(mode >= 0); + assert(mode < BC7CompressionMode::kNumModes); + if(gModeEstimate[mode] == -1.0 || est < gModeEstimate[mode]) { + gModeEstimate[mode] = est; + } } // Compress a single block. @@ -1825,12 +1872,31 @@ namespace BC7C } else { const float *w = GetErrorMetric(); - gModeEstimate[6] = 0.0001 + blockCluster.QuantizedError(Min, Max, 4, 0xFFF0F0F0, RGBAVector(w[0], w[1], w[2], w[3])); + const double err = 0.0001 + blockCluster.QuantizedError(Min, Max, 4, 0xFEFEFEFE, RGBAVector(w[0], w[1], w[2], w[3])); + UpdateErrorEstimate(6, err); + +#ifdef USE_PCA_FOR_SHAPE_ESTIMATION + if(statManager) { + double eigOne = blockCluster.GetPrincipalEigenvalue(); + double eigTwo = blockCluster.GetSecondEigenvalue(); + double error; + if(eigOne != 0.0) { + error = eigTwo / eigOne; + } + else { + error = 1.0; + } + + BlockStat s (kBlockStatString[eBlockStat_SingleShapeEstimate], error); + statManager->AddStat(blockIdx, s); + } +#endif } } // First we must figure out which shape to use. To do this, simply // see which shape has the smallest sum of minimum bounding spheres. + double estimates[2] = { -1.0, -1.0 }; double bestError[2] = { DBL_MAX, DBL_MAX }; int bestShapeIdx[2] = { -1, -1 }; RGBACluster bestClusters[2][3]; @@ -1841,8 +1907,38 @@ namespace BC7C PopulateTwoClustersForShape(blockCluster, i, clusters); double err = 0.0; + double errEstimate[2] = { -1.0, -1.0 }; for(int ci = 0; ci < 2; ci++) { - err += EstimateTwoClusterError(clusters[ci]); + double shapeEstimate[2] = { -1.0, -1.0 }; + err += EstimateTwoClusterError(clusters[ci], shapeEstimate); + + for(int ei = 0; ei < 2; ei++) { + if(shapeEstimate[ei] >= 0.0) { + if(errEstimate[ei] == -1.0) { + errEstimate[ei] = shapeEstimate[ei]; + } + else { + errEstimate[ei] += shapeEstimate[ei]; + } + } + } + } + +#ifdef USE_PCA_FOR_SHAPE_ESTIMATION + err /= 2.0; +#endif + + if(errEstimate[0] != -1.0) { + UpdateErrorEstimate(1, errEstimate[0]); + } + + if(errEstimate[1] != -1.0) { + UpdateErrorEstimate(3, errEstimate[1]); + } + + if(statManager && err < bestError[0]) { + BlockStat s = BlockStat(kBlockStatString[eBlockStat_TwoShapeEstimate], err); + statManager->AddStat(blockIdx, s); } // If it's small, we'll take it! @@ -1874,9 +1970,39 @@ namespace BC7C PopulateThreeClustersForShape(blockCluster, i, clusters); double err = 0.0; - for(int ci = 0; ci < 3; ci++) { - err += EstimateThreeClusterError(clusters[ci]); - } + double errEstimate[2] = { -1.0, -1.0 }; + for(int ci = 0; ci < 3; ci++) { + double shapeEstimate[2] = { -1.0, -1.0 }; + err += EstimateThreeClusterError(clusters[ci], shapeEstimate); + + for(int ei = 0; ei < 2; ei++) { + if(shapeEstimate[ei] >= 0.0) { + if(errEstimate[ei] == -1.0) { + errEstimate[ei] = shapeEstimate[ei]; + } + else { + errEstimate[ei] += shapeEstimate[ei]; + } + } + } + } + +#ifdef USE_PCA_FOR_SHAPE_ESTIMATION + err /= 3.0; +#endif + + if(errEstimate[0] != -1.0) { + UpdateErrorEstimate(0, errEstimate[0]); + } + + if(errEstimate[1] != -1.0) { + UpdateErrorEstimate(2, errEstimate[1]); + } + + if(statManager && err < bestError[1]) { + BlockStat s = BlockStat(kBlockStatString[eBlockStat_ThreeShapeEstimate], err); + statManager->AddStat(blockIdx, s); + } // If it's small, we'll take it! if(err < 1e-9) { diff --git a/BPTCEncoder/src/RGBAEndpoints.cpp b/BPTCEncoder/src/RGBAEndpoints.cpp index 7a52f6b..512af09 100755 --- a/BPTCEncoder/src/RGBAEndpoints.cpp +++ b/BPTCEncoder/src/RGBAEndpoints.cpp @@ -307,12 +307,52 @@ void RGBACluster::GetPrincipalAxis(RGBADir &axis) { } RGBAVector avg = m_Total / float(m_NumPoints); - ::GetPrincipalAxis(m_NumPoints, m_DataPoints, m_PrincipalAxis); + m_PowerMethodIterations = ::GetPrincipalAxis( + m_NumPoints, + m_DataPoints, + m_PrincipalAxis, + m_PrincipalEigenvalue, + &m_SecondEigenvalue + ); + m_PrincipalAxisCached = true; GetPrincipalAxis(axis); } +double RGBACluster::GetPrincipalEigenvalue() { + + if(!m_PrincipalAxisCached) { + RGBADir dummy; + GetPrincipalAxis(dummy); + } + + assert(m_PrincipalAxisCached); + return m_PrincipalEigenvalue; +} + +double RGBACluster::GetSecondEigenvalue() { + + if(!m_PrincipalAxisCached) { + RGBADir dummy; + GetPrincipalAxis(dummy); + } + + assert(m_PrincipalAxisCached); + return m_SecondEigenvalue; +} + +uint32 RGBACluster::GetPowerMethodIterations() { + + if(!m_PrincipalAxisCached) { + RGBADir dummy; + GetPrincipalAxis(dummy); + } + + assert(m_PrincipalAxisCached); + return m_PowerMethodIterations; +} + double RGBACluster::QuantizedError(const RGBAVector &p1, const RGBAVector &p2, uint8 nBuckets, uint32 bitMask, const RGBAVector &errorMetricVec, const int pbits[2], int *indices) const { // nBuckets should be a power of two. @@ -403,7 +443,58 @@ void ClampEndpoints(RGBAVector &p1, RGBAVector &p2) { clamp(p2.a, 0.0f, 255.0f); } -void GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis) { +static uint32 PowerIteration(const RGBAMatrix &mat, RGBADir &eigVec, double &eigVal) { + + int numIterations = 0; + const int kMaxNumIterations = 200; + + for(int nTries = 0; nTries < 3; nTries++) { + // !SPEED! Find eigenvectors by using the power method. This is good because the + // matrix is only 4x4, which allows us to use SIMD... + RGBAVector b = RGBAVector(rand()); + assert(b.Length() > 0); + b /= b.Length(); + + bool fixed = false; + numIterations = 0; + while(!fixed && ++numIterations < kMaxNumIterations) { + + RGBAVector newB = mat * b; + + // !HACK! If the principal eigenvector of the covariance matrix + // converges to zero, that means that the points lie equally + // spaced on a sphere in this space. In this (extremely rare) + // situation, just choose a point and use it as the principal + // direction. + const float newBlen = newB.Length(); + if(newBlen < 1e-10) { + eigVec = b; + eigVal = 0.0; + return numIterations; + } + + eigVal = newB.Length(); + newB /= eigVal; + + if(fabs(1.0f - (b * newB)) < 1e-5) + fixed = true; + + b = newB; + } + + eigVec = b; + if(numIterations < kMaxNumIterations) { + break; + } + } + + if(numIterations == kMaxNumIterations) { + eigVal = 0.0; + } + return numIterations; +} + +uint32 GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis, double &eigOne, double *eigTwo) { assert(nPts > 0); assert(nPts <= kMaxNumDataPoints); @@ -445,7 +536,7 @@ void GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis) { if(uptsIdx == 1) { axis.r = axis.g = axis.b = axis.a = 0.0f; - return; + return 0; } // Collinear? else { @@ -462,7 +553,7 @@ void GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis) { if(collinear) { axis = dir; - return; + return 0; } } @@ -481,39 +572,37 @@ void GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis) { covMatrix(j, i) = covMatrix(i, j); } } + + uint32 iters = PowerIteration(covMatrix, axis, eigOne); - // !SPEED! Find eigenvectors by using the power method. This is good because the - // matrix is only 4x4, which allows us to use SIMD... - RGBAVector b = toPtsMax; - assert(b.Length() > 0); - b /= b.Length(); + if(NULL != eigTwo) { + if(eigOne != 0.0) { + RGBAMatrix reduced = covMatrix - eigOne * RGBAMatrix( + axis.c[0] * axis.c[0], axis.c[0] * axis.c[1], axis.c[0] * axis.c[2], axis.c[0] * axis.c[3], + axis.c[1] * axis.c[0], axis.c[1] * axis.c[1], axis.c[1] * axis.c[2], axis.c[1] * axis.c[3], + axis.c[2] * axis.c[0], axis.c[2] * axis.c[1], axis.c[2] * axis.c[2], axis.c[2] * axis.c[3], + axis.c[3] * axis.c[0], axis.c[3] * axis.c[1], axis.c[3] * axis.c[2], axis.c[3] * axis.c[3] + ); + + bool allZero = true; + for(int i = 0; i < 16; i++) { + if(fabs(reduced[i]) > 0.0005) { + allZero = false; + } + } - bool fixed = false; - int infLoopPrevention = 0; - const int kMaxNumIterations = 200; - while(!fixed && ++infLoopPrevention < kMaxNumIterations) { + if(allZero) { + *eigTwo = 0.0; + } + else { + RGBADir dummyDir; + iters += PowerIteration(reduced, dummyDir, *eigTwo); + } + } + else { + *eigTwo = 0.0; + } + } - RGBAVector newB = covMatrix * b; - - // !HACK! If the principal eigenvector of the covariance matrix - // converges to zero, that means that the points lie equally - // spaced on a sphere in this space. In this (extremely rare) - // situation, just choose a point and use it as the principal - // direction. - const float newBlen = newB.Length(); - if(newBlen < 1e-10) { - axis = toPts[0]; - return; - } - - newB /= newB.Length(); - - if(fabs(1.0f - (b * newB)) < 1e-5) - fixed = true; - - b = newB; - } - - assert(infLoopPrevention < kMaxNumIterations); - axis = b; + return iters; } diff --git a/BPTCEncoder/src/RGBAEndpoints.h b/BPTCEncoder/src/RGBAEndpoints.h index 253b2a5..dcecd32 100755 --- a/BPTCEncoder/src/RGBAEndpoints.h +++ b/BPTCEncoder/src/RGBAEndpoints.h @@ -161,6 +161,18 @@ private: } public: + + RGBAMatrix( + float _m1, float _m2, float _m3, float _m4, + float _m5, float _m6, float _m7, float _m8, + float _m9, float _m10, float _m11, float _m12, + float _m13, float _m14, float _m15, float _m16 + ) : + m1(_m1), m2(_m2), m3(_m3), m4(_m4), + m5(_m5), m6(_m6), m7(_m7), m8(_m8), + m9(_m9), m10(_m10), m11(_m11), m12(_m12), + m13(_m13), m14(_m14), m15(_m15), m16(_m16) + { } RGBAMatrix() : m1(1.0f), m2(0.0f), m3(0.0f), m4(0.0f), @@ -294,7 +306,10 @@ public: m_PointBitString(c.m_PointBitString), m_Min(c.m_Min), m_Max(c.m_Max), - m_PrincipalAxisCached(false) + m_PrincipalAxisCached(c.m_PrincipalAxisCached), + m_SecondEigenvalue(c.m_SecondEigenvalue), + m_PowerMethodIterations(c.m_PowerMethodIterations), + m_PrincipalAxis(c.m_PrincipalAxis) { memcpy(this->m_DataPoints, c.m_DataPoints, m_NumPoints * sizeof(RGBAVector)); } @@ -327,6 +342,9 @@ public: double QuantizedError(const RGBAVector &p1, const RGBAVector &p2, uint8 nBuckets, uint32 bitMask, const RGBAVector &errorMetricVec, const int pbits[2] = NULL, int *indices = NULL) const; // Returns the principal axis for this point cluster. + double GetPrincipalEigenvalue(); + double GetSecondEigenvalue(); + uint32 GetPowerMethodIterations(); void GetPrincipalAxis(RGBADir &axis); bool AllSamePoint() const { return m_Max == m_Min; } @@ -345,11 +363,14 @@ private: RGBAVector m_Min, m_Max; int m_PointBitString; + double m_PrincipalEigenvalue; + double m_SecondEigenvalue; + uint32 m_PowerMethodIterations; RGBADir m_PrincipalAxis; bool m_PrincipalAxisCached; }; extern uint8 QuantizeChannel(const uint8 val, const uint8 mask, const int pBit = -1); -extern void GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis); +extern uint32 GetPrincipalAxis(int nPts, const RGBAVector *pts, RGBADir &axis, double &eigOne, double *eigTwo); #endif //__RGBA_ENDPOINTS_H__