- 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

@@ -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++) {
+	double errEstimate[2] = { -1.0, -1.0 };
-          err += EstimateThreeClusterError(clusters[ci]);
+	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) {

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
+	if(NULL != eigTwo) {
-	// matrix is only 4x4, which allows us to use SIMD...
+	  if(eigOne != 0.0) {
-	RGBAVector b = toPtsMax;
+	    RGBAMatrix reduced = covMatrix - eigOne * RGBAMatrix(
-	assert(b.Length() > 0);
+	      axis.c[0] * axis.c[0], axis.c[0] * axis.c[1], axis.c[0] * axis.c[2], axis.c[0] * axis.c[3], 
-	b /= b.Length();
+	      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;
+	    if(allZero) {
-	int infLoopPrevention = 0;
+	      *eigTwo = 0.0;
-	const int kMaxNumIterations = 200;
+	    }
-	while(!fixed && ++infLoopPrevention < kMaxNumIterations) {
+	    else {
+	      RGBADir dummyDir;
+	      iters += PowerIteration(reduced, dummyDir, *eigTwo);
+	    }
+	  }
+	  else {
+	    *eigTwo = 0.0;
+	  }
+        }
 
-		RGBAVector newB = covMatrix * b;
+	return iters;
-
-		// !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;
 }

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__