Add tool for aligning color space of video files

This class adds logic for aligning color space of a test video compared
to a reference video. If there is a color space mismatch, it typically
does not have much impact on human perception, but it has a big impact
on PSNR and SSIM calculations. For example, aligning a test run with VP8
improves PSNR and SSIM from:
Average PSNR: 29.142818, average SSIM: 0.946026
to:
Average PSNR: 38.146229, average SSIM: 0.965388.

The optiomal color transformation between the two videos were:
0.86 0.01 0.00 14.37
0.00 0.88 0.00 15.32
0.00 0.00 0.88 15.74
which is converting YUV full range to YUV limited range. There is
already a CL out for fixing this discrepancy here:
https://webrtc-review.googlesource.com/c/src/+/94543

After that, hopefully there is no color space mismatch when saving the
raw YUV values. It's good that the video quality tool is color space
agnostic anyway, and can compensate for differences when the test
video is obtained by e.g. filming a physical device screen.

Also, the linear least square logic will be used for compensating
geometric distorisions in a follow-up CL.

Bug: webrtc:9642
Change-Id: I499713960a0544d8e45c5d09886e68ec829b28a7
Reviewed-on: https://webrtc-review.googlesource.com/c/95950
Reviewed-by: Sami Kalliomäki <sakal@webrtc.org>
Reviewed-by: Patrik Höglund <phoglund@webrtc.org>
Commit-Queue: Magnus Jedvert <magjed@webrtc.org>
Cr-Commit-Position: refs/heads/master@{#25193}
diff --git a/rtc_tools/BUILD.gn b/rtc_tools/BUILD.gn
index c0716b4..3d9804c 100644
--- a/rtc_tools/BUILD.gn
+++ b/rtc_tools/BUILD.gn
@@ -86,6 +86,10 @@
 
 rtc_static_library("video_quality_analysis") {
   sources = [
+    "frame_analyzer/linear_least_squares.cc",
+    "frame_analyzer/linear_least_squares.h",
+    "frame_analyzer/video_color_aligner.cc",
+    "frame_analyzer/video_color_aligner.h",
     "frame_analyzer/video_quality_analysis.cc",
     "frame_analyzer/video_quality_analysis.h",
     "frame_analyzer/video_temporal_aligner.cc",
@@ -93,6 +97,7 @@
   ]
   deps = [
     ":video_file_reader",
+    "../api:array_view",
     "../api/video:video_frame_i420",
     "../common_video",
     "../rtc_base:checks",
@@ -333,7 +338,9 @@
     testonly = true
 
     sources = [
+      "frame_analyzer/linear_least_squares_unittest.cc",
       "frame_analyzer/reference_less_video_analysis_unittest.cc",
+      "frame_analyzer/video_color_aligner_unittest.cc",
       "frame_analyzer/video_quality_analysis_unittest.cc",
       "frame_analyzer/video_temporal_aligner_unittest.cc",
       "frame_editing/frame_editing_unittest.cc",
@@ -363,6 +370,7 @@
       "//test:test_support",
       "//testing/gtest",
       "//third_party/abseil-cpp/absl/memory",
+      "//third_party/libyuv",
     ]
 
     if (rtc_enable_protobuf) {
diff --git a/rtc_tools/frame_analyzer/linear_least_squares.cc b/rtc_tools/frame_analyzer/linear_least_squares.cc
new file mode 100644
index 0000000..706d212
--- /dev/null
+++ b/rtc_tools/frame_analyzer/linear_least_squares.cc
@@ -0,0 +1,200 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_tools/frame_analyzer/linear_least_squares.h"
+
+#include <numeric>
+#include <utility>
+
+#include "rtc_base/checks.h"
+#include "rtc_base/logging.h"
+
+namespace webrtc {
+namespace test {
+
+template <class T>
+using Matrix = std::valarray<std::valarray<T>>;
+
+namespace {
+
+template <typename R, typename T>
+R DotProduct(const std::valarray<T>& a, const std::valarray<T>& b) {
+  RTC_CHECK_EQ(a.size(), b.size());
+  return std::inner_product(std::begin(a), std::end(a), std::begin(b), R(0));
+}
+
+// Calculates a^T * b.
+template <typename R, typename T>
+Matrix<R> MatrixMultiply(const Matrix<T>& a, const Matrix<T>& b) {
+  Matrix<R> result(std::valarray<R>(a.size()), b.size());
+  for (size_t i = 0; i < a.size(); ++i) {
+    for (size_t j = 0; j < b.size(); ++j)
+      result[j][i] = DotProduct<R>(a[i], b[j]);
+  }
+
+  return result;
+}
+
+template <typename T>
+Matrix<T> Transpose(const Matrix<T>& matrix) {
+  if (matrix.size() == 0)
+    return Matrix<T>();
+  const size_t rows = matrix.size();
+  const size_t columns = matrix[0].size();
+  Matrix<T> result(std::valarray<T>(rows), columns);
+
+  for (size_t i = 0; i < rows; ++i) {
+    for (size_t j = 0; j < columns; ++j)
+      result[j][i] = matrix[i][j];
+  }
+
+  return result;
+}
+
+// Convert valarray from type T to type R.
+template <typename R, typename T>
+std::valarray<R> ConvertTo(const std::valarray<T>& v) {
+  std::valarray<R> result(v.size());
+  for (size_t i = 0; i < v.size(); ++i)
+    result[i] = static_cast<R>(v[i]);
+  return result;
+}
+
+// Convert valarray Matrix from type T to type R.
+template <typename R, typename T>
+Matrix<R> ConvertTo(const Matrix<T>& mat) {
+  Matrix<R> result(mat.size());
+  for (size_t i = 0; i < mat.size(); ++i)
+    result[i] = ConvertTo<R>(mat[i]);
+  return result;
+}
+
+// Convert from valarray Matrix back to the more conventional std::vector.
+template <typename T>
+std::vector<std::vector<T>> ToVectorMatrix(const Matrix<T>& m) {
+  std::vector<std::vector<T>> result;
+  for (const std::valarray<T>& v : m)
+    result.emplace_back(std::begin(v), std::end(v));
+  return result;
+}
+
+// Create a valarray Matrix from a conventional std::vector.
+template <typename T>
+Matrix<T> FromVectorMatrix(const std::vector<std::vector<T>>& mat) {
+  Matrix<T> result(mat.size());
+  for (size_t i = 0; i < mat.size(); ++i)
+    result[i] = std::valarray<T>(mat[i].data(), mat[i].size());
+  return result;
+}
+
+// Returns |matrix_to_invert|^-1 * |right_hand_matrix|. |matrix_to_invert| must
+// have square size.
+Matrix<double> GaussianElimination(Matrix<double> matrix_to_invert,
+                                   Matrix<double> right_hand_matrix) {
+  // |n| is the width/height of |matrix_to_invert|.
+  const size_t n = matrix_to_invert.size();
+  // Make sure |matrix_to_invert| has square size.
+  for (const std::valarray<double>& column : matrix_to_invert)
+    RTC_CHECK_EQ(n, column.size());
+  // Make sure |right_hand_matrix| has correct size.
+  for (const std::valarray<double>& column : right_hand_matrix)
+    RTC_CHECK_EQ(n, column.size());
+
+  // Transpose the matrices before and after so that we can perform Gaussian
+  // elimination on the columns instead of the rows, since that is easier with
+  // our representation.
+  matrix_to_invert = Transpose(matrix_to_invert);
+  right_hand_matrix = Transpose(right_hand_matrix);
+
+  // Loop over the diagonal of |matrix_to_invert| and perform column reduction.
+  // Column reduction is a sequence of elementary column operations that is
+  // performed on both |matrix_to_invert| and |right_hand_matrix| until
+  // |matrix_to_invert| has been transformed to the identity matrix.
+  for (size_t diagonal_index = 0; diagonal_index < n; ++diagonal_index) {
+    // Make sure the diagonal element has the highest absolute value by
+    // swapping columns if necessary.
+    for (size_t column = diagonal_index + 1; column < n; ++column) {
+      if (std::abs(matrix_to_invert[column][diagonal_index]) >
+          std::abs(matrix_to_invert[diagonal_index][diagonal_index])) {
+        std::swap(matrix_to_invert[column], matrix_to_invert[diagonal_index]);
+        std::swap(right_hand_matrix[column], right_hand_matrix[diagonal_index]);
+      }
+    }
+
+    // Reduce the diagonal element to be 1, by dividing the column with that
+    // value. If the diagonal element is 0, it means the system of equations has
+    // many solutions, and in that case we will return an arbitrary solution.
+    if (matrix_to_invert[diagonal_index][diagonal_index] == 0.0) {
+      RTC_LOG(LS_WARNING) << "Matrix is not invertible, ignoring.";
+      continue;
+    }
+    const double diagonal_element =
+        matrix_to_invert[diagonal_index][diagonal_index];
+    matrix_to_invert[diagonal_index] /= diagonal_element;
+    right_hand_matrix[diagonal_index] /= diagonal_element;
+
+    // Eliminate the other entries in row |diagonal_index| by making them zero.
+    for (size_t column = 0; column < n; ++column) {
+      if (column == diagonal_index)
+        continue;
+      const double row_element = matrix_to_invert[column][diagonal_index];
+      matrix_to_invert[column] -=
+          row_element * matrix_to_invert[diagonal_index];
+      right_hand_matrix[column] -=
+          row_element * right_hand_matrix[diagonal_index];
+    }
+  }
+
+  // Transpose the result before returning it, explained in comment above.
+  return Transpose(right_hand_matrix);
+}
+
+}  // namespace
+
+IncrementalLinearLeastSquares::IncrementalLinearLeastSquares() = default;
+IncrementalLinearLeastSquares::~IncrementalLinearLeastSquares() = default;
+
+void IncrementalLinearLeastSquares::AddObservations(
+    const std::vector<std::vector<uint8_t>>& x,
+    const std::vector<std::vector<uint8_t>>& y) {
+  if (x.empty() || y.empty())
+    return;
+  // Make sure all columns are the same size.
+  const size_t n = x[0].size();
+  for (const std::vector<uint8_t>& column : x)
+    RTC_CHECK_EQ(n, column.size());
+  for (const std::vector<uint8_t>& column : y)
+    RTC_CHECK_EQ(n, column.size());
+
+  // We will multiply the uint8_t values together, so we need to expand to a
+  // type that can safely store those values, i.e. uint16_t.
+  const Matrix<uint16_t> unpacked_x = ConvertTo<uint16_t>(FromVectorMatrix(x));
+  const Matrix<uint16_t> unpacked_y = ConvertTo<uint16_t>(FromVectorMatrix(y));
+
+  const Matrix<uint64_t> xx = MatrixMultiply<uint64_t>(unpacked_x, unpacked_x);
+  const Matrix<uint64_t> xy = MatrixMultiply<uint64_t>(unpacked_x, unpacked_y);
+  if (sum_xx && sum_xy) {
+    *sum_xx += xx;
+    *sum_xy += xy;
+  } else {
+    sum_xx = xx;
+    sum_xy = xy;
+  }
+}
+
+std::vector<std::vector<double>>
+IncrementalLinearLeastSquares::GetBestSolution() const {
+  RTC_CHECK(sum_xx && sum_xy) << "No observations have been added";
+  return ToVectorMatrix(GaussianElimination(ConvertTo<double>(*sum_xx),
+                                            ConvertTo<double>(*sum_xy)));
+}
+
+}  // namespace test
+}  // namespace webrtc
diff --git a/rtc_tools/frame_analyzer/linear_least_squares.h b/rtc_tools/frame_analyzer/linear_least_squares.h
new file mode 100644
index 0000000..1b07dc1
--- /dev/null
+++ b/rtc_tools/frame_analyzer/linear_least_squares.h
@@ -0,0 +1,53 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_
+#define RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_
+
+#include <valarray>
+#include <vector>
+#include "absl/types/optional.h"
+
+namespace webrtc {
+namespace test {
+
+// This class is used for finding a matrix b that roughly solves the equation:
+// y = x * b. This is generally impossible to do exactly, so the problem is
+// rephrased as finding the matrix b that minimizes the difference:
+// |y - x * b|^2. Calling multiple AddObservations() is equivalent to
+// concatenating the observation vectors and calling AddObservations() once. The
+// reason for doing it incrementally is that we can't store the raw YUV values
+// for a whole video file in memory at once. This class has a constant memory
+// footprint, regardless how may times AddObservations() is called.
+class IncrementalLinearLeastSquares {
+ public:
+  IncrementalLinearLeastSquares();
+  ~IncrementalLinearLeastSquares();
+
+  // Add a number of observations. The subvectors of x and y must have the same
+  // length.
+  void AddObservations(const std::vector<std::vector<uint8_t>>& x,
+                       const std::vector<std::vector<uint8_t>>& y);
+
+  // Calculate and return the best linear solution, given the observations so
+  // far.
+  std::vector<std::vector<double>> GetBestSolution() const;
+
+ private:
+  // Running sum of x^T * x.
+  absl::optional<std::valarray<std::valarray<uint64_t>>> sum_xx;
+  // Running sum of x^T * y.
+  absl::optional<std::valarray<std::valarray<uint64_t>>> sum_xy;
+};
+
+}  // namespace test
+}  // namespace webrtc
+
+#endif  // RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_
diff --git a/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc b/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc
new file mode 100644
index 0000000..b074aac
--- /dev/null
+++ b/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc
@@ -0,0 +1,91 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_tools/frame_analyzer/linear_least_squares.h"
+
+#include "test/gtest.h"
+
+namespace webrtc {
+namespace test {
+
+TEST(LinearLeastSquares, ScalarIdentityOneObservation) {
+  IncrementalLinearLeastSquares lls;
+  lls.AddObservations({{1}}, {{1}});
+  EXPECT_EQ(std::vector<std::vector<double>>({{1.0}}), lls.GetBestSolution());
+}
+
+TEST(LinearLeastSquares, ScalarIdentityTwoObservationsOneCall) {
+  IncrementalLinearLeastSquares lls;
+  lls.AddObservations({{1, 2}}, {{1, 2}});
+  EXPECT_EQ(std::vector<std::vector<double>>({{1.0}}), lls.GetBestSolution());
+}
+
+TEST(LinearLeastSquares, ScalarIdentityTwoObservationsTwoCalls) {
+  IncrementalLinearLeastSquares lls;
+  lls.AddObservations({{1}}, {{1}});
+  lls.AddObservations({{2}}, {{2}});
+  EXPECT_EQ(std::vector<std::vector<double>>({{1.0}}), lls.GetBestSolution());
+}
+
+TEST(LinearLeastSquares, MatrixIdentityOneObservation) {
+  IncrementalLinearLeastSquares lls;
+  lls.AddObservations({{1, 2}, {3, 4}}, {{1, 2}, {3, 4}});
+  EXPECT_EQ(std::vector<std::vector<double>>({{1.0, 0.0}, {0.0, 1.0}}),
+            lls.GetBestSolution());
+}
+
+TEST(LinearLeastSquares, MatrixManyObservations) {
+  IncrementalLinearLeastSquares lls;
+  // Test that we can find the solution of the overspecified equation system:
+  // [1, 2] [1, 3] = [5,  11]
+  // [3, 4] [2, 4]   [11, 25]
+  // [5, 6]          [17, 39]
+  lls.AddObservations({{1}, {2}}, {{5}, {11}});
+  lls.AddObservations({{3}, {4}}, {{11}, {25}});
+  lls.AddObservations({{5}, {6}}, {{17}, {39}});
+
+  const std::vector<std::vector<double>> result = lls.GetBestSolution();
+  // We allow some numerical flexibility here.
+  EXPECT_DOUBLE_EQ(1.0, result[0][0]);
+  EXPECT_DOUBLE_EQ(2.0, result[0][1]);
+  EXPECT_DOUBLE_EQ(3.0, result[1][0]);
+  EXPECT_DOUBLE_EQ(4.0, result[1][1]);
+}
+
+TEST(LinearLeastSquares, MatrixVectorOneObservation) {
+  IncrementalLinearLeastSquares lls;
+  // Test that we can find the solution of the overspecified equation system:
+  // [1, 2] [1] = [5]
+  // [3, 4] [2]   [11]
+  // [5, 6]       [17]
+  lls.AddObservations({{1, 3, 5}, {2, 4, 6}}, {{5, 11, 17}});
+
+  const std::vector<std::vector<double>> result = lls.GetBestSolution();
+  // We allow some numerical flexibility here.
+  EXPECT_DOUBLE_EQ(1.0, result[0][0]);
+  EXPECT_DOUBLE_EQ(2.0, result[0][1]);
+}
+
+TEST(LinearLeastSquares, LinearLeastSquaresNonPerfectSolution) {
+  IncrementalLinearLeastSquares lls;
+  // Test that we can find the non-perfect solution of the overspecified
+  // equation system:
+  // [1] [20] = [21]
+  // [2]        [39]
+  // [3]        [60]
+  // [2]        [41]
+  // [1]        [19]
+  lls.AddObservations({{1, 2, 3, 2, 1}}, {{21, 39, 60, 41, 19}});
+
+  EXPECT_DOUBLE_EQ(20.0, lls.GetBestSolution()[0][0]);
+}
+
+}  // namespace test
+}  // namespace webrtc
diff --git a/rtc_tools/frame_analyzer/video_color_aligner.cc b/rtc_tools/frame_analyzer/video_color_aligner.cc
new file mode 100644
index 0000000..7afb1e4
--- /dev/null
+++ b/rtc_tools/frame_analyzer/video_color_aligner.cc
@@ -0,0 +1,239 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_tools/frame_analyzer/video_color_aligner.h"
+
+#include <algorithm>
+#include <cmath>
+#include <deque>
+#include <limits>
+#include <vector>
+
+#include "api/array_view.h"
+#include "api/video/i420_buffer.h"
+#include "rtc_base/checks.h"
+#include "rtc_base/logging.h"
+#include "rtc_base/refcountedobject.h"
+#include "rtc_tools/frame_analyzer/linear_least_squares.h"
+#include "rtc_tools/frame_analyzer/video_quality_analysis.h"
+#include "third_party/libyuv/include/libyuv/compare.h"
+#include "third_party/libyuv/include/libyuv/planar_functions.h"
+#include "third_party/libyuv/include/libyuv/scale.h"
+
+namespace webrtc {
+namespace test {
+
+namespace {
+
+// Helper function for AdjustColors(). This functions calculates a single output
+// row for y with the given color coefficients. The u/v channels are assumed to
+// be subsampled by a factor of 2, which is the case of I420.
+void CalculateYChannel(rtc::ArrayView<const uint8_t> y_data,
+                       rtc::ArrayView<const uint8_t> u_data,
+                       rtc::ArrayView<const uint8_t> v_data,
+                       const std::array<float, 4>& coeff,
+                       rtc::ArrayView<uint8_t> output) {
+  RTC_CHECK_EQ(y_data.size(), output.size());
+  // Each u/v element represents two y elements. Make sure we have enough to
+  // cover the Y values.
+  RTC_CHECK_GE(u_data.size() * 2, y_data.size());
+  RTC_CHECK_GE(v_data.size() * 2, y_data.size());
+
+  // Do two pixels at a time since u/v are subsampled.
+  for (size_t i = 0; i * 2 < y_data.size() - 1; ++i) {
+    const float uv_contribution =
+        coeff[1] * u_data[i] + coeff[2] * v_data[i] + coeff[3];
+
+    const float val0 = coeff[0] * y_data[i * 2 + 0] + uv_contribution;
+    const float val1 = coeff[0] * y_data[i * 2 + 1] + uv_contribution;
+
+    // Clamp result to a byte.
+    output[i * 2 + 0] = static_cast<uint8_t>(
+        std::round(std::max(0.0f, std::min(val0, 255.0f))));
+    output[i * 2 + 1] = static_cast<uint8_t>(
+        std::round(std::max(0.0f, std::min(val1, 255.0f))));
+  }
+
+  // Handle the last pixel for odd widths.
+  if (y_data.size() % 2 == 1) {
+    const float val = coeff[0] * y_data[y_data.size() - 1] +
+                      coeff[1] * u_data[(y_data.size() - 1) / 2] +
+                      coeff[2] * v_data[(y_data.size() - 1) / 2] + coeff[3];
+    output[y_data.size() - 1] =
+        static_cast<uint8_t>(std::round(std::max(0.0f, std::min(val, 255.0f))));
+  }
+}
+
+// Helper function for AdjustColors(). This functions calculates a single output
+// row for either u or v, with the given color coefficients. Y, U, and V are
+// assumed to be the same size, i.e. no subsampling.
+void CalculateUVChannel(rtc::ArrayView<const uint8_t> y_data,
+                        rtc::ArrayView<const uint8_t> u_data,
+                        rtc::ArrayView<const uint8_t> v_data,
+                        const std::array<float, 4>& coeff,
+                        rtc::ArrayView<uint8_t> output) {
+  RTC_CHECK_EQ(y_data.size(), u_data.size());
+  RTC_CHECK_EQ(y_data.size(), v_data.size());
+  RTC_CHECK_EQ(y_data.size(), output.size());
+
+  for (size_t x = 0; x < y_data.size(); ++x) {
+    const float val = coeff[0] * y_data[x] + coeff[1] * u_data[x] +
+                      coeff[2] * v_data[x] + coeff[3];
+    // Clamp result to a byte.
+    output[x] =
+        static_cast<uint8_t>(std::round(std::max(0.0f, std::min(val, 255.0f))));
+  }
+}
+
+// Convert a frame to four vectors consisting of [y, u, v, 1].
+std::vector<std::vector<uint8_t>> FlattenYuvData(
+    const rtc::scoped_refptr<I420BufferInterface>& frame) {
+  std::vector<std::vector<uint8_t>> result(
+      4, std::vector<uint8_t>(frame->ChromaWidth() * frame->ChromaHeight()));
+
+  // Downscale the Y plane so that all YUV planes are the same size.
+  libyuv::ScalePlane(frame->DataY(), frame->StrideY(), frame->width(),
+                     frame->height(), result[0].data(), frame->ChromaWidth(),
+                     frame->ChromaWidth(), frame->ChromaHeight(),
+                     libyuv::kFilterBox);
+
+  libyuv::CopyPlane(frame->DataU(), frame->StrideU(), result[1].data(),
+                    frame->ChromaWidth(), frame->ChromaWidth(),
+                    frame->ChromaHeight());
+
+  libyuv::CopyPlane(frame->DataV(), frame->StrideV(), result[2].data(),
+                    frame->ChromaWidth(), frame->ChromaWidth(),
+                    frame->ChromaHeight());
+
+  std::fill(result[3].begin(), result[3].end(), 1u);
+
+  return result;
+}
+
+ColorTransformationMatrix VectorToColorMatrix(
+    const std::vector<std::vector<double>>& v) {
+  ColorTransformationMatrix color_transformation;
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 4; ++j)
+      color_transformation[i][j] = v[i][j];
+  }
+  return color_transformation;
+}
+
+}  // namespace
+
+ColorTransformationMatrix CalculateColorTransformationMatrix(
+    const rtc::scoped_refptr<I420BufferInterface>& reference_frame,
+    const rtc::scoped_refptr<I420BufferInterface>& test_frame) {
+  IncrementalLinearLeastSquares incremental_lls;
+  incremental_lls.AddObservations(FlattenYuvData(test_frame),
+                                  FlattenYuvData(reference_frame));
+  return VectorToColorMatrix(incremental_lls.GetBestSolution());
+}
+
+ColorTransformationMatrix CalculateColorTransformationMatrix(
+    const rtc::scoped_refptr<Video>& reference_video,
+    const rtc::scoped_refptr<Video>& test_video) {
+  RTC_CHECK_GE(reference_video->number_of_frames(),
+               test_video->number_of_frames());
+
+  IncrementalLinearLeastSquares incremental_lls;
+  for (size_t i = 0; i < test_video->number_of_frames(); ++i) {
+    incremental_lls.AddObservations(
+        FlattenYuvData(test_video->GetFrame(i)),
+        FlattenYuvData(reference_video->GetFrame(i)));
+  }
+
+  return VectorToColorMatrix(incremental_lls.GetBestSolution());
+}
+
+rtc::scoped_refptr<Video> AdjustColors(
+    const ColorTransformationMatrix& color_transformation,
+    const rtc::scoped_refptr<Video>& video) {
+  class ColorAdjustedVideo : public rtc::RefCountedObject<Video> {
+   public:
+    ColorAdjustedVideo(const ColorTransformationMatrix& color_transformation,
+                       const rtc::scoped_refptr<Video>& video)
+        : color_transformation_(color_transformation), video_(video) {}
+
+    int width() const override { return video_->width(); }
+    int height() const override { return video_->height(); }
+    size_t number_of_frames() const override {
+      return video_->number_of_frames();
+    }
+
+    rtc::scoped_refptr<I420BufferInterface> GetFrame(
+        size_t index) const override {
+      return AdjustColors(color_transformation_, video_->GetFrame(index));
+    }
+
+   private:
+    const ColorTransformationMatrix color_transformation_;
+    const rtc::scoped_refptr<Video> video_;
+  };
+
+  return new ColorAdjustedVideo(color_transformation, video);
+}
+
+rtc::scoped_refptr<I420BufferInterface> AdjustColors(
+    const ColorTransformationMatrix& color_matrix,
+    const rtc::scoped_refptr<I420BufferInterface>& frame) {
+  // Allocate I420 buffer that will hold the color adjusted frame.
+  rtc::scoped_refptr<I420Buffer> adjusted_frame =
+      I420Buffer::Create(frame->width(), frame->height());
+
+  // Create a downscaled Y plane with the same size as the U/V planes to
+  // simplify converting the U/V planes.
+  std::vector<uint8_t> downscaled_y_plane(frame->ChromaWidth() *
+                                          frame->ChromaHeight());
+  libyuv::ScalePlane(frame->DataY(), frame->StrideY(), frame->width(),
+                     frame->height(), downscaled_y_plane.data(),
+                     frame->ChromaWidth(), frame->ChromaWidth(),
+                     frame->ChromaHeight(), libyuv::kFilterBox);
+
+  // Fill in the adjusted data row by row.
+  for (int y = 0; y < frame->height(); ++y) {
+    const int half_y = y / 2;
+    rtc::ArrayView<const uint8_t> y_row(frame->DataY() + frame->StrideY() * y,
+                                        frame->width());
+    rtc::ArrayView<const uint8_t> u_row(
+        frame->DataU() + frame->StrideU() * half_y, frame->ChromaWidth());
+    rtc::ArrayView<const uint8_t> v_row(
+        frame->DataV() + frame->StrideV() * half_y, frame->ChromaWidth());
+    rtc::ArrayView<uint8_t> output_y_row(
+        adjusted_frame->MutableDataY() + adjusted_frame->StrideY() * y,
+        frame->width());
+
+    CalculateYChannel(y_row, u_row, v_row, color_matrix[0], output_y_row);
+
+    // Chroma channels only exist every second row for I420.
+    if (y % 2 == 0) {
+      rtc::ArrayView<const uint8_t> downscaled_y_row(
+          downscaled_y_plane.data() + frame->ChromaWidth() * half_y,
+          frame->ChromaWidth());
+      rtc::ArrayView<uint8_t> output_u_row(
+          adjusted_frame->MutableDataU() + adjusted_frame->StrideU() * half_y,
+          frame->ChromaWidth());
+      rtc::ArrayView<uint8_t> output_v_row(
+          adjusted_frame->MutableDataV() + adjusted_frame->StrideV() * half_y,
+          frame->ChromaWidth());
+
+      CalculateUVChannel(downscaled_y_row, u_row, v_row, color_matrix[1],
+                         output_u_row);
+      CalculateUVChannel(downscaled_y_row, u_row, v_row, color_matrix[2],
+                         output_v_row);
+    }
+  }
+
+  return adjusted_frame;
+}
+
+}  // namespace test
+}  // namespace webrtc
diff --git a/rtc_tools/frame_analyzer/video_color_aligner.h b/rtc_tools/frame_analyzer/video_color_aligner.h
new file mode 100644
index 0000000..997cbf6
--- /dev/null
+++ b/rtc_tools/frame_analyzer/video_color_aligner.h
@@ -0,0 +1,49 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef RTC_TOOLS_FRAME_ANALYZER_VIDEO_COLOR_ALIGNER_H_
+#define RTC_TOOLS_FRAME_ANALYZER_VIDEO_COLOR_ALIGNER_H_
+
+#include <array>
+
+#include "rtc_tools/video_file_reader.h"
+
+namespace webrtc {
+namespace test {
+
+// Represents a linear color transformation from [y, u, v] to [y', u', v']
+// through the equation: [y', u', v'] = [y, u, v, 1] * matrix.
+using ColorTransformationMatrix = std::array<std::array<float, 4>, 3>;
+
+// Calculate the optimal color transformation that should be applied to the test
+// video to match as closely as possible to the reference video.
+ColorTransformationMatrix CalculateColorTransformationMatrix(
+    const rtc::scoped_refptr<Video>& reference_video,
+    const rtc::scoped_refptr<Video>& test_video);
+
+// Calculate color transformation for a single I420 frame.
+ColorTransformationMatrix CalculateColorTransformationMatrix(
+    const rtc::scoped_refptr<I420BufferInterface>& reference_frame,
+    const rtc::scoped_refptr<I420BufferInterface>& test_frame);
+
+// Apply a color transformation to a video.
+rtc::scoped_refptr<Video> AdjustColors(
+    const ColorTransformationMatrix& color_matrix,
+    const rtc::scoped_refptr<Video>& video);
+
+// Apply a color transformation to a single I420 frame.
+rtc::scoped_refptr<I420BufferInterface> AdjustColors(
+    const ColorTransformationMatrix& color_matrix,
+    const rtc::scoped_refptr<I420BufferInterface>& frame);
+
+}  // namespace test
+}  // namespace webrtc
+
+#endif  // RTC_TOOLS_FRAME_ANALYZER_VIDEO_COLOR_ALIGNER_H_
diff --git a/rtc_tools/frame_analyzer/video_color_aligner_unittest.cc b/rtc_tools/frame_analyzer/video_color_aligner_unittest.cc
new file mode 100644
index 0000000..520d6f4
--- /dev/null
+++ b/rtc_tools/frame_analyzer/video_color_aligner_unittest.cc
@@ -0,0 +1,174 @@
+/*
+ *  Copyright (c) 2018 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_tools/frame_analyzer/video_color_aligner.h"
+
+#include "rtc_tools/frame_analyzer/video_quality_analysis.h"
+#include "rtc_tools/video_file_reader.h"
+#include "test/gtest.h"
+#include "test/testsupport/fileutils.h"
+#include "third_party/libyuv/include/libyuv/scale.h"
+
+namespace webrtc {
+namespace test {
+
+namespace {
+
+const ColorTransformationMatrix kIdentityColorMatrix = {
+    {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}}};
+
+void ExpectNear(const ColorTransformationMatrix& expected,
+                const ColorTransformationMatrix& actual) {
+  // The scaling factor on y/u/v should be pretty precise.
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j)
+      EXPECT_NEAR(expected[i][j], actual[i][j], /* abs_error= */ 1.0e-3)
+          << "at element i: " << i << ", j: " << j;
+  }
+  // The offset can be less precise since the range is [0, 255].
+  for (int i = 0; i < 3; ++i)
+    EXPECT_NEAR(expected[i][3], actual[i][3], /* abs_error= */ 0.1)
+        << "at element i: " << i;
+}
+
+}  // namespace
+
+class VideoColorAlignerTest : public ::testing::Test {
+ protected:
+  void SetUp() {
+    reference_video_ =
+        OpenYuvFile(ResourcePath("foreman_128x96", "yuv"), 128, 96);
+    ASSERT_TRUE(reference_video_);
+  }
+
+  rtc::scoped_refptr<Video> reference_video_;
+};
+
+TEST_F(VideoColorAlignerTest, AdjustColorsFrameIdentity) {
+  const rtc::scoped_refptr<I420BufferInterface> test_frame =
+      reference_video_->GetFrame(0);
+
+  // Assume perfect match, i.e. ssim == 1.
+  EXPECT_EQ(1.0,
+            Ssim(test_frame, AdjustColors(kIdentityColorMatrix, test_frame)));
+}
+
+TEST_F(VideoColorAlignerTest, AdjustColorsFrame1x1) {
+  const ColorTransformationMatrix color_matrix = {
+      {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}};
+
+  const uint8_t data_y[] = {2};
+  const uint8_t data_u[] = {6};
+  const uint8_t data_v[] = {7};
+  const rtc::scoped_refptr<I420BufferInterface> i420_buffer = I420Buffer::Copy(
+      /* width= */ 1, /* height= */ 1, data_y, /* stride_y= */ 1, data_u,
+      /* stride_u= */ 1, data_v, /* stride_v= */ 1);
+
+  const rtc::scoped_refptr<I420BufferInterface> adjusted_buffer =
+      AdjustColors(color_matrix, i420_buffer);
+
+  EXPECT_EQ(2 * 1 + 6 * 2 + 7 * 3 + 4, adjusted_buffer->DataY()[0]);
+  EXPECT_EQ(2 * 5 + 6 * 6 + 7 * 7 + 8, adjusted_buffer->DataU()[0]);
+  EXPECT_EQ(2 * 9 + 6 * 10 + 7 * 11 + 12, adjusted_buffer->DataV()[0]);
+}
+
+TEST_F(VideoColorAlignerTest, AdjustColorsFrame1x1Negative) {
+  const ColorTransformationMatrix color_matrix = {
+      {{-1, 0, 0, 255}, {0, -1, 0, 255}, {0, 0, -1, 255}}};
+
+  const uint8_t data_y[] = {2};
+  const uint8_t data_u[] = {6};
+  const uint8_t data_v[] = {7};
+  const rtc::scoped_refptr<I420BufferInterface> i420_buffer = I420Buffer::Copy(
+      /* width= */ 1, /* height= */ 1, data_y, /* stride_y= */ 1, data_u,
+      /* stride_u= */ 1, data_v, /* stride_v= */ 1);
+
+  const rtc::scoped_refptr<I420BufferInterface> adjusted_buffer =
+      AdjustColors(color_matrix, i420_buffer);
+
+  EXPECT_EQ(255 - 2, adjusted_buffer->DataY()[0]);
+  EXPECT_EQ(255 - 6, adjusted_buffer->DataU()[0]);
+  EXPECT_EQ(255 - 7, adjusted_buffer->DataV()[0]);
+}
+
+TEST_F(VideoColorAlignerTest, AdjustColorsFrame2x2) {
+  const ColorTransformationMatrix color_matrix = {
+      {{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}}};
+
+  const uint8_t data_y[] = {0, 1, 3, 4};
+  const uint8_t data_u[] = {6};
+  const uint8_t data_v[] = {7};
+  const rtc::scoped_refptr<I420BufferInterface> i420_buffer = I420Buffer::Copy(
+      /* width= */ 2, /* height= */ 2, data_y, /* stride_y= */ 2, data_u,
+      /* stride_u= */ 1, data_v, /* stride_v= */ 1);
+
+  const rtc::scoped_refptr<I420BufferInterface> adjusted_buffer =
+      AdjustColors(color_matrix, i420_buffer);
+
+  EXPECT_EQ(0 * 1 + 6 * 2 + 7 * 3 + 4, adjusted_buffer->DataY()[0]);
+  EXPECT_EQ(1 * 1 + 6 * 2 + 7 * 3 + 4, adjusted_buffer->DataY()[1]);
+  EXPECT_EQ(3 * 1 + 6 * 2 + 7 * 3 + 4, adjusted_buffer->DataY()[2]);
+  EXPECT_EQ(4 * 1 + 6 * 2 + 7 * 3 + 4, adjusted_buffer->DataY()[3]);
+
+  EXPECT_EQ(2 * 5 + 6 * 6 + 7 * 7 + 8, adjusted_buffer->DataU()[0]);
+  EXPECT_EQ(2 * 9 + 6 * 10 + 7 * 11 + 12, adjusted_buffer->DataV()[0]);
+}
+
+TEST_F(VideoColorAlignerTest, CalculateColorTransformationMatrixIdentity) {
+  EXPECT_EQ(kIdentityColorMatrix, CalculateColorTransformationMatrix(
+                                      reference_video_, reference_video_));
+}
+
+TEST_F(VideoColorAlignerTest, CalculateColorTransformationMatrixOffset) {
+  const uint8_t small_data_y[] = {0, 1, 2,  3,  4,  5,  6,  7,
+                                  8, 9, 10, 11, 12, 13, 14, 15};
+  const uint8_t small_data_u[] = {15, 13, 17, 29};
+  const uint8_t small_data_v[] = {3, 200, 170, 29};
+  const rtc::scoped_refptr<I420BufferInterface> small_i420_buffer =
+      I420Buffer::Copy(
+          /* width= */ 4, /* height= */ 4, small_data_y, /* stride_y= */ 4,
+          small_data_u, /* stride_u= */ 2, small_data_v, /* stride_v= */ 2);
+
+  uint8_t big_data_y[16];
+  uint8_t big_data_u[4];
+  uint8_t big_data_v[4];
+  // Create another I420 frame where all values are 10 bigger.
+  for (int i = 0; i < 16; ++i)
+    big_data_y[i] = small_data_y[i] + 10;
+  for (int i = 0; i < 4; ++i)
+    big_data_u[i] = small_data_u[i] + 10;
+  for (int i = 0; i < 4; ++i)
+    big_data_v[i] = small_data_v[i] + 10;
+
+  const rtc::scoped_refptr<I420BufferInterface> big_i420_buffer =
+      I420Buffer::Copy(
+          /* width= */ 4, /* height= */ 4, big_data_y, /* stride_y= */ 4,
+          big_data_u, /* stride_u= */ 2, big_data_v, /* stride_v= */ 2);
+
+  const ColorTransformationMatrix color_matrix =
+      CalculateColorTransformationMatrix(big_i420_buffer, small_i420_buffer);
+
+  ExpectNear({{{1, 0, 0, 10}, {0, 1, 0, 10}, {0, 0, 1, 10}}}, color_matrix);
+}
+
+TEST_F(VideoColorAlignerTest, CalculateColorTransformationMatrix) {
+  // Arbitrary color transformation matrix.
+  const ColorTransformationMatrix org_color_matrix = {
+      {{0.8, 0.05, 0.04, -4}, {-0.2, 0.7, 0.1, 10}, {0.1, 0.2, 0.4, 20}}};
+
+  const ColorTransformationMatrix result_color_matrix =
+      CalculateColorTransformationMatrix(
+          AdjustColors(org_color_matrix, reference_video_), reference_video_);
+
+  ExpectNear(org_color_matrix, result_color_matrix);
+}
+
+}  // namespace test
+}  // namespace webrtc