New Matrix3F class added in gfx.

ui/gfx/matrix3_f.cc

Index: ui/gfx/matrix3_f.cc
diff --git a/ui/gfx/matrix3_f.cc b/ui/gfx/matrix3_f.cc
new file mode 100644
index 0000000000000000000000000000000000000000..72088dfb35b46732c8b13e610eac0e3181e950e2
--- /dev/null
+++ b/ui/gfx/matrix3_f.cc
@@ -0,0 +1,301 @@
+// 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"
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+namespace {
+
+// This is only to make accessing indices self-explanatory.
+enum MatrixCoordinates {
+  kM00,

[danakj, 2013/01/31 22:38:30]

enum names should be CAPS without the k prefix.

[motek., 2013/01/31 23:29:56]

Is this a recent change to style guide? The only other enum declared in a cc
file in this module (in skbitmap_operation.cc) uses the k-camel convention. 

[danakj, 2013/01/31 23:32:56]

There was a recent thread on chromium-dev a month or so back about it. There is
a bit of inconsistency in existing code, but the consensus which is in the style
guide now was:

"Though the Google C++ Style Guide now says to use kConstantNaming for enums,
Chromium was written using MACRO_STYLE naming.  Continue to use this style for
consistency."

+  kM01,
+  kM02,
+  kM10,
+  kM11,
+  kM12,
+  kM20,
+  kM21,
+  kM22,
+  kMEnd
+};
+
+template<typename T>
+double Determinant3x3(T data[kMEnd]) {
+  // 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[kM00]) * (
+          static_cast<double>(data[kM11]) * data[kM22] -
+          static_cast<double>(data[kM12]) * data[kM21]) +
+      static_cast<double>(data[kM01]) * (
+          static_cast<double>(data[kM12]) * data[kM20] -
+          static_cast<double>(data[kM10]) * data[kM22]) +
+      static_cast<double>(data[kM02]) * (
+          static_cast<double>(data[kM10]) * data[kM21] -
+          static_cast<double>(data[kM11]) * data[kM20]);
+}
+
+int MatrixToArrayCoords(int i, int j) {

[danakj, 2013/01/31 22:38:30]

make this a private member function in the header so you can use it from get()
and set()

[motek., 2013/01/31 23:29:56]

Done.

+  DCHECK(i >= 0 && i < 3);
+  DCHECK(j >= 0 && j < 3);
+  return i * 3 + j;
+}
+
+}  // namespace
+
+namespace gfx {
+
+Matrix3F::Matrix3F() {
+}
+
+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 {
+  return data_[MatrixToArrayCoords(i, j)];

[danakj, 2013/01/31 22:38:30]

move this to the header

[motek., 2013/01/31 23:29:56]

Done.

+}
+
+void Matrix3F::set(int i, int j, ScalarType v) {
+  data_[MatrixToArrayCoords(i, j)] = v;

[danakj, 2013/01/31 22:38:30]

move this to the header

[motek., 2013/01/31 23:29:56]

Done.

+}
+
+void Matrix3F::set(ScalarType v) {

[danakj, 2013/01/31 22:38:30]

set the whole matrix to a single value? is this really something that will be
done commonly? i've never seen this on a matrix before.

[motek., 2013/01/31 23:29:56]

Honestly? I had no use for it either. It has been requested by the other
asvitkine@ in the review (used from Zeros() and Ones()). 
Should I get rid of it? Sure, its gone. 

[Alexei Svitkine (slow), 2013/01/31 23:34:11]

Deferring to Dana here, I just suggested it because it's more concise from the
call site than specifying 9 constants. But if that's an uncommon use, then kill
it.

+  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 {

[danakj, 2013/01/31 22:38:30]

move to header. get_column()

[motek., 2013/01/31 23:29:56]

Done.

+  DCHECK(i >= 0 && i < 3);
+  return VectorType(data_[i], data_[3 + i], data_[6 + i]);

[Alexei Svitkine (slow), 2013/01/31 22:17:18]

Use MatrixToArrayCoords() here.

[motek., 2013/01/31 22:28:59]

Done.

+}
+
+void Matrix3F::SetColumn(int i, const VectorType& c) {

[danakj, 2013/01/31 22:38:30]

move to header

[motek., 2013/01/31 23:29:56]

Done.

+  DCHECK(i >= 0 && i < 3);
+  data_[i] = c.x();

[Alexei Svitkine (slow), 2013/01/31 22:17:18]

Use MatrixToArrayCoords() here.

[motek., 2013/01/31 22:28:59]

Done.

+  data_[3 + i] = c.y();
+  data_[6 + i] = c.z();
+}
+
+Matrix3F Matrix3F::Inverse() const {
+  Matrix3F inverse = Matrix3F::Zeros();
+  double determinant = Determinant3x3(data_);
+  if (std::numeric_limits<ScalarType>::epsilon() > std::abs(determinant))
+    return inverse;  // Singular matrix. Return Zeros().
+
+  inverse.set(
+      (data_[kM11] * data_[kM22] - data_[kM12] * data_[kM21]) / determinant,
+      (data_[kM02] * data_[kM21] - data_[kM01] * data_[kM22]) / determinant,
+      (data_[kM01] * data_[kM12] - data_[kM02] * data_[kM11]) / determinant,
+      (data_[kM12] * data_[kM20] - data_[kM10] * data_[kM22]) / determinant,
+      (data_[kM00] * data_[kM22] - data_[kM02] * data_[kM20]) / determinant,
+      (data_[kM02] * data_[kM10] - data_[kM00] * data_[kM12]) / determinant,
+      (data_[kM10] * data_[kM21] - data_[kM11] * data_[kM20]) / determinant,
+      (data_[kM01] * data_[kM20] - data_[kM00] * data_[kM21]) / determinant,
+      (data_[kM00] * data_[kM11] - data_[kM01] * data_[kM10]) / determinant);
+  return inverse;
+}
+
+Matrix3F::ScalarType Matrix3F::Determinant() const {
+  return ScalarType(Determinant3x3(data_));

[danakj, 2013/01/31 22:38:30]

static cast<ScalarType>()

[motek., 2013/01/31 23:29:56]

Done.

+}
+
+Matrix3F::ScalarType Matrix3F::Trace() const {

[danakj, 2013/01/31 22:38:30]

move to header

[motek., 2013/01/31 23:29:56]

Done.

+  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);

[danakj, 2013/01/31 22:38:30]

if this is the only use of set(float), i'd prefer set(0,0,0,0,0,0,0,0,0) here
and in Ones()

[motek., 2013/01/31 23:29:56]

Done.

+  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);

[danakj, 2013/01/31 22:38:30]

use 1.f or 1.0f for float literals. same for 0.0, throughout this file.

[motek., 2013/01/31 23:29:56]

Done.

+  return matrix;
+}
+
+// static
+Matrix3F Matrix3F::FromOuterProduct(const VectorType& a,
+                                    const VectorType& 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;
+}
+
+}  // namespace gfx