New Matrix3F class added in gfx.

ui/gfx/matrix3_f.cc

+// 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>
+
+#include "base/logging.h"
+
+namespace {
+
+// This is only to make accessing indices self-explanatory.
+enum MatrixCoordinates {
+  kM00,
+  kM01,
+  kM02,
+  kM10,
+  kM11,
+  kM12,
+  kM20,
+  kM21,
+  kM22,
+  kMEnd
+};
+
+template<typename T>
+double Determinant3x3(T data[kMEnd]) {

[Alexei Svitkine (slow), 2013/01/31 20:02:38]

Add a comment explaining the purpose (same as you reply to me that you need a
higher precision version.)

[motek., 2013/01/31 22:00:24]

Done.

+  return (double)data[kM00] * ((double)data[kM11] * data[kM22] -

[Alexei Svitkine (slow), 2013/01/31 20:02:38]

Use C++ casts, rather than C-style casts.

[motek., 2013/01/31 22:00:24]

Done.

+                               (double)data[kM12] * data[kM21]) +
+      data[kM01] * ((double)data[kM12] * data[kM20] -
+                    (double)data[kM10] * data[kM22]) +
+      data[kM02] * ((double)data[kM10] * data[kM21] -
+                    (double)data[kM11] * data[kM20]);
+}
+
+}  // namespace
+
+namespace gfx {
+
+Matrix3F::Matrix3F() {
+}
+
+Matrix3F::Matrix3F(const VectorType& a, const VectorType& bt) {
+  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());
+}
+
+Matrix3F::Matrix3F(const Matrix3F& rhs) {
+  memcpy(data_, rhs.data_, sizeof(data_));
+}
+
+Matrix3F::~Matrix3F() {
+}
+
+bool Matrix3F::IsEqual(const Matrix3F& rhs) const {
+  return 0 == memcmp(data_, rhs.data_, sizeof(data_));
+}
+
+bool Matrix3F::IsNear(const Matrix3F& rhs, ScalarType precision) const {
+  DCHECK(precision >= 0);
+  for (int i = 0; i < kMEnd; ++i) {
+    if (std::abs(data_[i] - rhs.data_[i]) > precision)
+      return false;
+  }
+  return true;
+}
+
+Matrix3F::ScalarType Matrix3F::get(int i, int j) const {
+  DCHECK(i >= 0 && i < 3);
+  DCHECK(j >= 0 && j < 3);
+  return data_[i * 3 + j];

[Alexei Svitkine (slow), 2013/01/31 20:02:38]

Make a function in the anon namespace that gives you the index, and stick the
dchecks there. Then, use that whenever you index into the array via i,j.

e.g.

int GetElementIndex(int i, int j) {
  DCHECK(i >= 0 && i < 3);
  DCHECK(j >= 0 && j < 3);
  return i * 3 + j;
}

[motek., 2013/01/31 22:00:24]

Done.

+}
+
+void Matrix3F::set(int i, int j, ScalarType v) {
+  DCHECK(i >= 0 && i < 3);
+  DCHECK(j >= 0 && j < 3);
+  data_[i * 3 + j] = v;
+}
+
+void Matrix3F::set(ScalarType v) {
+  for (int i = 0; i < kMEnd; i++)
+    data_[i] = v;
+}
+
+void Matrix3F::set(ScalarType m00, ScalarType m01, ScalarType m02,
+                   ScalarType m10, ScalarType m11, ScalarType m12,
+                   ScalarType m20, ScalarType m21, ScalarType 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;
+}
+
+Matrix3F::VectorType Matrix3F::GetColumn(int i) const {
+  DCHECK(i >= 0 && i < 3);
+  return VectorType(data_[i], data_[3 + i], data_[6 + i]);
+}
+
+void Matrix3F::SetColumn(int i, const VectorType& c) {
+  DCHECK(i >= 0 && i < 3);
+  data_[i] = c.x();
+  data_[3 + i] = c.y();
+  data_[6 + i] = c.z();
+}
+
+Matrix3F Matrix3F::Inverse() const {
+  Matrix3F inverse = Matrix3F::Zeros();
+  double det = Determinant3x3(data_);

[Alexei Svitkine (slow), 2013/01/31 20:02:38]

Nit: det -> determinant

[motek., 2013/01/31 22:00:24]

Done.

+  if (std::numeric_limits<ScalarType>::epsilon() > std::abs(det))
+    return inverse;  // Singular matrix. Return Zeros().
+
+  inverse.set((data_[kM11] * data_[kM22] - data_[kM12] * data_[kM21]) / det,
+              (data_[kM02] * data_[kM21] - data_[kM01] * data_[kM22]) / det,
+              (data_[kM01] * data_[kM12] - data_[kM02] * data_[kM11]) / det,
+              (data_[kM12] * data_[kM20] - data_[kM10] * data_[kM22]) / det,
+              (data_[kM00] * data_[kM22] - data_[kM02] * data_[kM20]) / det,
+              (data_[kM02] * data_[kM10] - data_[kM00] * data_[kM12]) / det,
+              (data_[kM10] * data_[kM21] - data_[kM11] * data_[kM20]) / det,
+              (data_[kM01] * data_[kM20] - data_[kM00] * data_[kM21]) / det,
+              (data_[kM00] * data_[kM11] - data_[kM01] * data_[kM10]) / det);
+  return inverse;
+}
+
+Matrix3F::ScalarType Matrix3F::Determinant() const {
+  return ScalarType(Determinant3x3(data_));
+}
+
+Matrix3F::ScalarType Matrix3F::Trace() const {
+  return data_[kM00] + data_[kM11] + data_[kM22];
+}
+
+Matrix3F::VectorType Matrix3F::SolveEigenproblem(
+    Matrix3F* eigenvectors) const {
+  // The matrix must be symmetric.
+  const ScalarType epsilon = std::numeric_limits<ScalarType>::epsilon();
+  if (std::abs(data_[kM01] - data_[kM10]) > epsilon ||
+      std::abs(data_[kM02] - data_[kM02]) > epsilon ||
+      std::abs(data_[kM12] - data_[kM21]) > epsilon) {
+    NOTREACHED();
+    return Matrix3F::VectorType();
+  }
+
+  ScalarType eigenvalues[3];
+  ScalarType p =
+      data_[kM01] * data_[kM01] +
+      data_[kM02] * data_[kM02] +
+      data_[kM12] * data_[kM12];
+
+  bool diagonal = std::abs(p) < epsilon;
+  if (diagonal) {
+    eigenvalues[0] = data_[kM00];
+    eigenvalues[1] = data_[kM11];
+    eigenvalues[2] = data_[kM22];
+  } else {
+    ScalarType q = Trace() / 3.0;
+    p = (data_[kM00] - q) * (data_[kM00] - q) +
+        (data_[kM11] - q) * (data_[kM11] - q) +
+        (data_[kM22] - q) * (data_[kM22] - q) +
+        2 * p;
+    p = sqrt(p / 6);
+
+    // The computation below puts B as (A - qI) / p, where A is *this.
+    Matrix3F matrix_B(*this);
+    matrix_B.data_[kM00] -= q;
+    matrix_B.data_[kM11] -= q;
+    matrix_B.data_[kM22] -= q;
+    for (int i = 0; i < kMEnd; ++i)
+      matrix_B.data_[i] /= p;
+
+    ScalarType half_det_b = matrix_B.Determinant() / 2.0;
+    // half_det_b should be in <-1, 1>, but beware of rounding error.
+    ScalarType phi = 0.0;
+    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 * cos(phi);
+    eigenvalues[2] = q + 2 * p * 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.0);
+  } 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) {
+      ScalarType l = eigenvalues[i];
+      // B = A - l * I
+      Matrix3F matrix_B(*this);
+      matrix_B.data_[kM00] -= l;
+      matrix_B.data_[kM11] -= l;
+      matrix_B.data_[kM22] -= l;
+      VectorType e1 = CrossProduct(matrix_B.GetColumn(0),
+                                   matrix_B.GetColumn(1));
+      VectorType e2 = CrossProduct(matrix_B.GetColumn(1),
+                                   matrix_B.GetColumn(2));
+      VectorType e3 = CrossProduct(matrix_B.GetColumn(2),
+                                   matrix_B.GetColumn(0));
+
+      // e1, e2 and e3 should point in the same direction.
+      if (DotProduct(e1, e2) < 0) {
+        e2.set_x(-e2.x());
+        e2.set_y(-e2.y());
+        e2.set_z(-e2.z());
+      }
+
+      if (DotProduct(e1, e3) < 0) {
+        e3.set_x(-e3.x());
+        e3.set_y(-e3.y());
+        e3.set_z(-e3.z());
+      }
+
+      VectorType eigvec = e1 + e2 + e3;
+      // Normalize.
+      eigvec.Scale(1.0 / eigvec.Length());
+      eigenvectors->SetColumn(i, eigvec);
+    }
+  }
+
+  return Matrix3F::VectorType(eigenvalues[0],
+                              eigenvalues[1],
+                              eigenvalues[2]);
+}
+
+// static
+Matrix3F Matrix3F::Zeros() {
+  Matrix3F matrix;
+  matrix.set(0.0);
+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::Ones() {
+  Matrix3F matrix;
+  matrix.set(1.0);
+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::Identity() {
+  Matrix3F matrix;
+  matrix.set(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
+  return matrix;
+}
+
+}  // namespace gfx