Added Matrix3F class to gfx.

Will need it to implement color decomposition for thumbnails.


BUG=155269

Review URL: https://chromiumcodereview.appspot.com/12096069

git-svn-id: svn://svn.chromium.org/chrome/trunk/src@180427 0039d316-1c4b-4281-b951-d872f2087c98

3 files changed, 493 insertions, 0 deletions

diff --git a/ui/gfx/matrix3_f.cc b/ui/gfx/matrix3_f.cc
new file mode 100644
index 0000000..501c5dd
--- /dev/null
+++ b/ui/gfx/matrix3_f.cc
@@ -0,0 +1,237 @@
+// Copyright (c) 2013 The Chromium Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include "ui/gfx/matrix3_f.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+namespace {
+
+// This is only to make accessing indices self-explanatory.
+enum MatrixCoordinates {
+  M00,
+  M01,
+  M02,
+  M10,
+  M11,
+  M12,
+  M20,
+  M21,
+  M22,
+  M_END
+};
+
+template<typename T>
+double Determinant3x3(T data[M_END]) {
+  // This routine is separated from the Matrix3F::Determinant because in
+  // computing inverse we do want higher precision afforded by the explicit
+  // use of 'double'.
+  return
+      static_cast<double>(data[M00]) * (
+          static_cast<double>(data[M11]) * data[M22] -
+          static_cast<double>(data[M12]) * data[M21]) +
+      static_cast<double>(data[M01]) * (
+          static_cast<double>(data[M12]) * data[M20] -
+          static_cast<double>(data[M10]) * data[M22]) +
+      static_cast<double>(data[M02]) * (
+          static_cast<double>(data[M10]) * data[M21] -
+          static_cast<double>(data[M11]) * data[M20]);
+}
+
+}  // namespace
+
+namespace gfx {
+
+Matrix3F::Matrix3F() {
+}
+
+Matrix3F::~Matrix3F() {
+}
+
+// static
+Matrix3F Matrix3F::Zeros() {
+  Matrix3F matrix;
+  matrix.set(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::Ones() {
+  Matrix3F matrix;
+  matrix.set(1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f);
+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::Identity() {
+  Matrix3F matrix;
+  matrix.set(1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f);
+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::FromOuterProduct(const Vector3dF& a, const Vector3dF& bt) {
+  Matrix3F matrix;
+  matrix.set(a.x() * bt.x(), a.x() * bt.y(), a.x() * bt.z(),
+             a.y() * bt.x(), a.y() * bt.y(), a.y() * bt.z(),
+             a.z() * bt.x(), a.z() * bt.y(), a.z() * bt.z());
+  return matrix;
+}
+
+bool Matrix3F::IsEqual(const Matrix3F& rhs) const {
+  return 0 == memcmp(data_, rhs.data_, sizeof(data_));
+}
+
+bool Matrix3F::IsNear(const Matrix3F& rhs, float precision) const {
+  DCHECK(precision >= 0);
+  for (int i = 0; i < M_END; ++i) {
+    if (std::abs(data_[i] - rhs.data_[i]) > precision)
+      return false;
+  }
+  return true;
+}
+
+Matrix3F Matrix3F::Inverse() const {
+  Matrix3F inverse = Matrix3F::Zeros();
+  double determinant = Determinant3x3(data_);
+  if (std::numeric_limits<float>::epsilon() > std::abs(determinant))
+    return inverse;  // Singular matrix. Return Zeros().
+
+  inverse.set(
+      (data_[M11] * data_[M22] - data_[M12] * data_[M21]) / determinant,
+      (data_[M02] * data_[M21] - data_[M01] * data_[M22]) / determinant,
+      (data_[M01] * data_[M12] - data_[M02] * data_[M11]) / determinant,
+      (data_[M12] * data_[M20] - data_[M10] * data_[M22]) / determinant,
+      (data_[M00] * data_[M22] - data_[M02] * data_[M20]) / determinant,
+      (data_[M02] * data_[M10] - data_[M00] * data_[M12]) / determinant,
+      (data_[M10] * data_[M21] - data_[M11] * data_[M20]) / determinant,
+      (data_[M01] * data_[M20] - data_[M00] * data_[M21]) / determinant,
+      (data_[M00] * data_[M11] - data_[M01] * data_[M10]) / determinant);
+  return inverse;
+}
+
+float Matrix3F::Determinant() const {
+  return static_cast<float>(Determinant3x3(data_));
+}
+
+Vector3dF Matrix3F::SolveEigenproblem(Matrix3F* eigenvectors) const {
+  // The matrix must be symmetric.
+  const float epsilon = std::numeric_limits<float>::epsilon();
+  if (std::abs(data_[M01] - data_[M10]) > epsilon ||
+      std::abs(data_[M02] - data_[M02]) > epsilon ||
+      std::abs(data_[M12] - data_[M21]) > epsilon) {
+    NOTREACHED();
+    return Vector3dF();
+  }
+
+  float eigenvalues[3];
+  float p =
+      data_[M01] * data_[M01] +
+      data_[M02] * data_[M02] +
+      data_[M12] * data_[M12];
+
+  bool diagonal = std::abs(p) < epsilon;
+  if (diagonal) {
+    eigenvalues[0] = data_[M00];
+    eigenvalues[1] = data_[M11];
+    eigenvalues[2] = data_[M22];
+  } else {
+    float q = Trace() / 3.0f;
+    p = (data_[M00] - q) * (data_[M00] - q) +
+        (data_[M11] - q) * (data_[M11] - q) +
+        (data_[M22] - q) * (data_[M22] - q) +
+        2 * p;
+    p = std::sqrt(p / 6);
+
+    // The computation below puts B as (A - qI) / p, where A is *this.
+    Matrix3F matrix_b(*this);
+    matrix_b.data_[M00] -= q;
+    matrix_b.data_[M11] -= q;
+    matrix_b.data_[M22] -= q;
+    for (int i = 0; i < M_END; ++i)
+      matrix_b.data_[i] /= p;
+
+    double half_det_b = Determinant3x3(matrix_b.data_) / 2.0;
+    // half_det_b should be in <-1, 1>, but beware of rounding error.
+    double phi = 0.0f;
+    if (half_det_b <= -1.0)
+      phi = M_PI / 3;
+    else if (half_det_b < 1.0)
+      phi = acos(half_det_b) / 3;
+
+    eigenvalues[0] = q + 2 * p * static_cast<float>(cos(phi));
+    eigenvalues[2] = q + 2 * p *
+        static_cast<float>(cos(phi + 2.0 * M_PI / 3.0));
+    eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
+  }
+
+  // Put eigenvalues in the descending order.
+  int indices[3] = {0, 1, 2};
+  if (eigenvalues[2] > eigenvalues[1]) {
+    std::swap(eigenvalues[2], eigenvalues[1]);
+    std::swap(indices[2], indices[1]);
+  }
+
+  if (eigenvalues[1] > eigenvalues[0]) {
+    std::swap(eigenvalues[1], eigenvalues[0]);
+    std::swap(indices[1], indices[0]);
+  }
+
+  if (eigenvalues[2] > eigenvalues[1]) {
+    std::swap(eigenvalues[2], eigenvalues[1]);
+    std::swap(indices[2], indices[1]);
+  }
+
+  if (eigenvectors != NULL && diagonal) {
+    // Eigenvectors are e-vectors, just need to be sorted accordingly.
+    *eigenvectors = Zeros();
+    for (int i = 0; i < 3; ++i)
+      eigenvectors->set(indices[i], i, 1.0f);
+  } else if (eigenvectors != NULL) {
+    // Consult the following for a detailed discussion:
+    // Joachim Kopp
+    // Numerical diagonalization of hermitian 3x3 matrices
+    // arXiv.org preprint: physics/0610206
+    // Int. J. Mod. Phys. C19 (2008) 523-548
+
+    // TODO(motek): expand to handle correctly negative and multiple
+    // eigenvalues.
+    for (int i = 0; i < 3; ++i) {
+      float l = eigenvalues[i];
+      // B = A - l * I
+      Matrix3F matrix_b(*this);
+      matrix_b.data_[M00] -= l;
+      matrix_b.data_[M11] -= l;
+      matrix_b.data_[M22] -= l;
+      Vector3dF e1 = CrossProduct(matrix_b.get_column(0),
+                                  matrix_b.get_column(1));
+      Vector3dF e2 = CrossProduct(matrix_b.get_column(1),
+                                  matrix_b.get_column(2));
+      Vector3dF e3 = CrossProduct(matrix_b.get_column(2),
+                                  matrix_b.get_column(0));
+
+      // e1, e2 and e3 should point in the same direction.
+      if (DotProduct(e1, e2) < 0)
+        e2 = -e2;
+
+      if (DotProduct(e1, e3) < 0)
+        e3 = -e3;
+
+      Vector3dF eigvec = e1 + e2 + e3;
+      // Normalize.
+      eigvec.Scale(1.0f / eigvec.Length());
+      eigenvectors->set_column(i, eigvec);
+    }
+  }
+
+  return Vector3dF(eigenvalues[0], eigenvalues[1], eigenvalues[2]);
+}
+
+}  // namespace gfx
diff --git a/ui/gfx/matrix3_f.h b/ui/gfx/matrix3_f.h
new file mode 100644
index 0000000..83f9cd9
--- /dev/null
+++ b/ui/gfx/matrix3_f.h
@@ -0,0 +1,108 @@
+// Copyright (c) 2013 The Chromium Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#ifndef UI_GFX_MATRIX3_F_H_
+#define UI_GFX_MATRIX3_F_H_
+
+#include "base/logging.h"
+#include "ui/gfx/vector3d_f.h"
+
+namespace gfx {
+
+class UI_EXPORT Matrix3F {
+ public:
+  ~Matrix3F();
+
+  static Matrix3F Zeros();
+  static Matrix3F Ones();
+  static Matrix3F Identity();
+  static Matrix3F FromOuterProduct(const Vector3dF& a, const Vector3dF& bt);
+
+  bool IsEqual(const Matrix3F& rhs) const;
+
+  // Element-wise comparison with given precision.
+  bool IsNear(const Matrix3F& rhs, float precision) const;
+
+  float get(int i, int j) const {
+    return data_[MatrixToArrayCoords(i, j)];
+  }
+
+  void set(int i, int j, float v) {
+    data_[MatrixToArrayCoords(i, j)] = v;
+  }
+
+  void set(float m00, float m01, float m02,
+           float m10, float m11, float m12,
+           float m20, float m21, float m22) {
+    data_[0] = m00;
+    data_[1] = m01;
+    data_[2] = m02;
+    data_[3] = m10;
+    data_[4] = m11;
+    data_[5] = m12;
+    data_[6] = m20;
+    data_[7] = m21;
+    data_[8] = m22;
+  }
+
+  Vector3dF get_column(int i) const {
+    return Vector3dF(
+      data_[MatrixToArrayCoords(0, i)],
+      data_[MatrixToArrayCoords(1, i)],
+      data_[MatrixToArrayCoords(2, i)]);
+  }
+
+  void set_column(int i, const Vector3dF& c) {
+    data_[MatrixToArrayCoords(0, i)] = c.x();
+    data_[MatrixToArrayCoords(1, i)] = c.y();
+    data_[MatrixToArrayCoords(2, i)] = c.z();
+  }
+
+  // Returns an inverse of this if the matrix is non-singular, zero (== Zero())
+  // otherwise.
+  Matrix3F Inverse() const;
+
+  // Value of the determinant of the matrix.
+  float Determinant() const;
+
+  // Trace (sum of diagonal elements) of the matrix.
+  float Trace() const {
+    return data_[MatrixToArrayCoords(0, 0)] +
+        data_[MatrixToArrayCoords(1, 1)] +
+        data_[MatrixToArrayCoords(2, 2)];
+  }
+
+  // Compute eigenvalues and (optionally) normalized eigenvectors of
+  // a positive defnite matrix *this. Eigenvectors are computed only if
+  // non-null |eigenvectors| matrix is passed. If it is NULL, the routine
+  // will not attempt to compute eigenvectors but will still return eigenvalues
+  // if they can be computed.
+  // If eigenvalues cannot be computed (the matrix does not meet constraints)
+  // the 0-vector is returned. Note that to retrieve eigenvalues, the matrix
+  // only needs to be symmetric while eigenvectors require it to be
+  // positive-definite. Passing a non-positive definite matrix will result in
+  // NaNs in vectors which cannot be computed.
+  // Eigenvectors are placed as column in |eigenvectors| in order corresponding
+  // to eigenvalues.
+  Vector3dF SolveEigenproblem(Matrix3F* eigenvectors) const;
+
+ private:
+  Matrix3F();  // Uninitialized default.
+
+  static int MatrixToArrayCoords(int i, int j) {
+    DCHECK(i >= 0 && i < 3);
+    DCHECK(j >= 0 && j < 3);
+    return i * 3 + j;
+  }
+
+  float data_[9];
+};
+
+inline bool operator==(const Matrix3F& lhs, const Matrix3F& rhs) {
+  return lhs.IsEqual(rhs);
+}
+
+}  // namespace gfx
+
+#endif  // UI_GFX_MATRIX3_F_H_
diff --git a/ui/gfx/matrix3_unittest.cc b/ui/gfx/matrix3_unittest.cc
new file mode 100644
index 0000000..04c656f
--- /dev/null
+++ b/ui/gfx/matrix3_unittest.cc
@@ -0,0 +1,148 @@
+// Copyright (c) 2013 The Chromium Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include <cmath>
+#include <limits>
+
+#include "base/basictypes.h"
+#include "testing/gtest/include/gtest/gtest.h"
+#include "ui/gfx/matrix3_f.h"
+
+namespace gfx {
+namespace {
+
+TEST(Matrix3fTest, Constructors) {
+  Matrix3F zeros = Matrix3F::Zeros();
+  Matrix3F ones = Matrix3F::Ones();
+  Matrix3F identity = Matrix3F::Identity();
+
+  Matrix3F product_ones = Matrix3F::FromOuterProduct(
+      Vector3dF(1.0f, 1.0f, 1.0f), Vector3dF(1.0f, 1.0f, 1.0f));
+  Matrix3F product_zeros = Matrix3F::FromOuterProduct(
+      Vector3dF(1.0f, 1.0f, 1.0f), Vector3dF(0.0f, 0.0f, 0.0f));
+  EXPECT_EQ(ones, product_ones);
+  EXPECT_EQ(zeros, product_zeros);
+
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j)
+      EXPECT_EQ(i == j ? 1.0f : 0.0f, identity.get(i, j));
+  }
+}
+
+TEST(Matrix3fTest, DataAccess) {
+  Matrix3F matrix = Matrix3F::Ones();
+  Matrix3F identity = Matrix3F::Identity();
+
+  EXPECT_EQ(Vector3dF(0.0f, 1.0f, 0.0f), identity.get_column(1));
+  matrix.set(0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f);
+  EXPECT_EQ(Vector3dF(2.0f, 5.0f, 8.0f), matrix.get_column(2));
+  matrix.set_column(0, Vector3dF(0.1f, 0.2f, 0.3f));
+  EXPECT_EQ(Vector3dF(0.1f, 0.2f, 0.3f), matrix.get_column(0));
+
+  EXPECT_EQ(0.1f, matrix.get(0, 0));
+  EXPECT_EQ(5.0f, matrix.get(1, 2));
+}
+
+TEST(Matrix3fTest, Determinant) {
+  EXPECT_EQ(1.0f, Matrix3F::Identity().Determinant());
+  EXPECT_EQ(0.0f, Matrix3F::Zeros().Determinant());
+  EXPECT_EQ(0.0f, Matrix3F::Ones().Determinant());
+
+  // Now for something non-trivial...
+  Matrix3F matrix = Matrix3F::Zeros();
+  matrix.set(0, 5, 6, 8, 7, 0, 1, 9, 0);
+  EXPECT_EQ(390.0f, matrix.Determinant());
+  matrix.set(2, 0, 3 * matrix.get(0, 0));
+  matrix.set(2, 1, 3 * matrix.get(0, 1));
+  matrix.set(2, 2, 3 * matrix.get(0, 2));
+  EXPECT_EQ(0, matrix.Determinant());
+
+  matrix.set(0.57f,  0.205f,  0.942f,
+             0.314f,  0.845f,  0.826f,
+             0.131f,  0.025f,  0.962f);
+  EXPECT_NEAR(0.3149f, matrix.Determinant(), 0.0001f);
+}
+
+TEST(Matrix3fTest, Inverse) {
+  Matrix3F identity = Matrix3F::Identity();
+  Matrix3F inv_identity = identity.Inverse();
+  EXPECT_EQ(identity, inv_identity);
+
+  Matrix3F singular = Matrix3F::Zeros();
+  singular.set(1.0f, 3.0f, 4.0f,
+               2.0f, 11.0f, 5.0f,
+               0.5f, 1.5f, 2.0f);
+  EXPECT_EQ(0, singular.Determinant());
+  EXPECT_EQ(Matrix3F::Zeros(), singular.Inverse());
+
+  Matrix3F regular = Matrix3F::Zeros();
+  regular.set(0.57f,  0.205f,  0.942f,
+              0.314f,  0.845f,  0.826f,
+              0.131f,  0.025f,  0.962f);
+  Matrix3F inv_regular = regular.Inverse();
+  regular.set(2.51540616f, -0.55138018f, -1.98968043f,
+              -0.61552266f,  1.34920184f, -0.55573636f,
+              -0.32653861f,  0.04002158f,  1.32488726f);
+  EXPECT_TRUE(regular.IsNear(inv_regular, 0.00001f));
+}
+
+TEST(Matrix3fTest, EigenvectorsIdentity) {
+  // This block tests the trivial case of eigenvalues of the identity matrix.
+  Matrix3F identity = Matrix3F::Identity();
+  Vector3dF eigenvals = identity.SolveEigenproblem(NULL);
+  EXPECT_EQ(Vector3dF(1.0f, 1.0f, 1.0f), eigenvals);
+}
+
+TEST(Matrix3fTest, EigenvectorsDiagonal)  {
+  // This block tests the another trivial case of eigenvalues of a diagonal
+  // matrix. Here we expect values to be sorted.
+  Matrix3F matrix = Matrix3F::Zeros();
+  matrix.set(0, 0, 1.0f);
+  matrix.set(1, 1, -2.5f);
+  matrix.set(2, 2, 3.14f);
+  Matrix3F eigenvectors = Matrix3F::Zeros();
+  Vector3dF eigenvals = matrix.SolveEigenproblem(&eigenvectors);
+  EXPECT_EQ(Vector3dF(3.14f, 1.0f, -2.5f), eigenvals);
+
+  EXPECT_EQ(Vector3dF(0.0f, 0.0f, 1.0f), eigenvectors.get_column(0));
+  EXPECT_EQ(Vector3dF(1.0f, 0.0f, 0.0f), eigenvectors.get_column(1));
+  EXPECT_EQ(Vector3dF(0.0f, 1.0f, 0.0f), eigenvectors.get_column(2));
+}
+
+TEST(Matrix3fTest, EigenvectorsNiceNotPositive)  {
+  // This block tests computation of eigenvectors of a matrix where nice
+  // round values are expected.
+  Matrix3F matrix = Matrix3F::Zeros();
+  // This is not a positive-definite matrix but eigenvalues and the first
+  // eigenvector should nonetheless be computed correctly.
+  matrix.set(3, 2, 4, 2, 0, 2, 4, 2, 3);
+  Matrix3F eigenvectors = Matrix3F::Zeros();
+  Vector3dF eigenvals = matrix.SolveEigenproblem(&eigenvectors);
+  EXPECT_EQ(Vector3dF(8.0f, -1.0f, -1.0f), eigenvals);
+
+  Vector3dF expected_principal(0.66666667f, 0.33333333f, 0.66666667f);
+  EXPECT_NEAR(0.0f,
+              (expected_principal - eigenvectors.get_column(0)).Length(),
+              0.000001f);
+}
+
+TEST(Matrix3fTest, EigenvectorsPositiveDefinite) {
+  // This block tests computation of eigenvectors of a matrix where output
+  // is not as nice as above, but it actually meets the definition.
+  Matrix3F matrix = Matrix3F::Zeros();
+  Matrix3F eigenvectors = Matrix3F::Zeros();
+  Matrix3F expected_eigenvectors = Matrix3F::Zeros();
+  matrix.set(1, -1,  2, -1,  4,  5, 2,  5,  0);
+  Vector3dF eigenvals =  matrix.SolveEigenproblem(&eigenvectors);
+  Vector3dF expected_eigv(7.3996266f, 1.91197255f, -4.31159915f);
+  expected_eigv -= eigenvals;
+  EXPECT_NEAR(0, expected_eigv.LengthSquared(), 0.00001f);
+  expected_eigenvectors.set(0.04926317f, -0.92135662f, -0.38558414f,
+                            0.82134249f, 0.25703273f, -0.50924521f,
+                            0.56830419f, -0.2916096f, 0.76941158f);
+  EXPECT_TRUE(expected_eigenvectors.IsNear(eigenvectors, 0.00001f));
+}
+
+}
+}